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ABSTRACT 


l 

Using  linear  systems  theory  as  a  framework,  the  solution  for  the  acoustic  field 
present  in  a  range -independent  acoustic  channel  excited  by  a  complex -weighted, 
planar  array  of  point  sources  with  an  arbitrary  input  electrical  signal  is  derived. 
The  ocean  medium  is  characterized  by  a  transfer  function,  obtainable  as  the  solution 
to  the  Helmholtz  wave  equation.  The  transfer  function  for  an  isospeed,  three-layer 
waveguide  is  derived.  The  unbounded  homogeneous  medium  equations  are  derived 
as  a  special  case  of  the  waveguide  problem.  The  problem  of  interference  due  to  the 
presence  of  a  pressure -release  surface  is  also  derived  as  a  special  case.  The  linear 
systems  approach  lends  itself  to  a  modular  computer  implementation,  in  which 
different  ocean  medium  models  are  represented  by  subroutines  implementing  their 
transfer  functions.  The  equations  for  a  range -independent  medium  are  implemented 
as  a  group  of  subprograms.  Results  are  presented  for  the  special  cases  of  a 
homogeneous  medium  and  the  surface  reflection  problem,  which  can  be  checked 
against  known,  easily  interpreted  analytical  solutions.  Finally,  an  example  of 
waveform  prediction  for  the  isospeed,  three-layer  waveguide  is  presented. 
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I.  INTRODUCTION 


Model -based  or  matched -field  signal  processing  refers  to  the  use  of  the 
knowledge  of  the  physics  of  a  problem  for  the  construction  of  mathematical  models 
and  its  application  to  suitable  signal  processing  algorithms.  Its  application  to 
underwater  acoustics  has  been  the  subject  of  many  papers  in  the  signal  processing 
literature,  as  in  Ziomek  and  Blount  [Ref.  1],  Baggeroer,  Kuperman  and 
Schmidt  [Ref.  2]  and  references  therein.  As  the  interest  in  this  field  grows  and  new 
algorithms  are  developed,  so  will  the  need  for  simulation  tools  to  verify  the  behavior 
of  those  algorithms.  Such  tools  should  be  able  to  generate  signals  with  spatial  and 
temporal  structures  like  those  found  in  real  environments. 

The  waveform  prediction  problem  has  been  studied  by  Officer  [Ref.  3:pp. 
101  -1 10,130]  for  the  specific  cases  of  pulse  reflection  from  a  boundary  and 
transmission  in  shallow  water.  DiNapoli  and  Deavenport  [Ref.4:pp.  131-134] 
describe  briefly  a  method  employing  the  Fast -Field -Program  (FFF)  suitable  to 
range -independent  problems.  Our  purpose  is  to  derive  a  generic  solution  for  the 
range-independent,  time-invariant,  deterministic  acoustic  channel  excited  bv  an 
array  of  point  sources  transmitting  an  arbitrary  signal.  In  our  derivation  we  will 
use  a  linear  systems  approach  to  the  acoustic  problem,  as  described  in 
Ziomek  [Ref.  5],  whose  notation  we  follow  closely.  In  particular,  we  will  apply  our 
results  to  the  isospeed  Pekeris  waveguide.  As  we  will  show,  the  linear  systems 
approach  lends  itself  to  a  modular  computer  implementation,  where  different  ocean 
media  are  represented  by  subroutines  or  functions  implementing  their  transfer 
function. 
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In  Section  II  we  present  the  basic  input -output  relationships  for  a 
space -variant,  time-independent  linear  system,  which  are  the  basis  of  our 
derivations.  Next,  we  show  how  those  equations  can  be  used  to  represent  an 
acoustic  channel.  The  particular  case  of  a  range -independent  ocean  is  then  studied, 
when  we  derive  the  output  waveform  equation  and  describe  the  process  of  derivation 
of  the  medium  transfer  function.  Finally,  we  apply  the  above  results  to  the  isospeed 
Pekeris  waveguide.  In  Section  III  we  discuss  computer  implementation  issues  and 
present  some  results  of  a  computer  implementation  of  the  Pekeris  waveguide 
equations. 


II.  ANALYSIS 


A.  THE  ACOUSTIC  CHANNEL  AS  A  LINEAR  SYSTEM 

In  this  section  we  will  present  some  results  from  the  theory  of  space -time  linear 
systems  [Ref.  5],  establish  the  notation  to  be  used,  and  derive  some  intermediate 
results  to  be  used  in  subsequent  sections. 

1.  Space-Variant,  Time -Invariant  Linear  Systems 

a.  Impulse  Response 

Linear,  time-invariant,  space-variant  filters,  as  shown  in  Fig.  1,  are 
represented  by  linear  partial  differential  equations  whose  coefficients  are  time 
independent.  In  linear  systems  terminology,  such  systems  are  characterized  by  their 
time-invariant,  space -dependent  impulse  response  A( r,r0;  r).  The  output  -f(tj)  is 
given  by  [Ref.  .ri;p.  o] 


rV-  r)  = 


a 


g{  t-T.r-fo)  h(  r,ru;r)  dr  dr. 


1) 


0 

J 


A  physical  interpretation  of  h ,  r  and  r0  can  be  obtained  by  letting  the 
input  be  an  impulse  applied  at  time  t  =  to  and  at  position  r  =  rs,  that  is, 
<?(t,r)  =  6(  t-t0, r-rs)  =  6{t-to)  5(r-rs).  Substituting  this  expression  for  g(  t.r)  into 
Eq.  (1)  and  using  the  sifting  property  of  the  impulse  function  we  obtain,  for  the 
output, 

<p{t,r)  =  h{  t-t0  ,  r-rs  ;  r  ),  (2) 

XX  XX 

r  r0 


from  which  we  can  define  h(r,r0;r)  as  the  response  of  the  system  at  time  t  and 
position  r  due  to  the  application  of  an  impulse  r  seconds  ago,  that  is,  at  the  instant 
to  =  f-r,  when  the  point  source  is  positioned  at  rs  =  r-r0.  Figure  2  illustrates  the 
geometry  of  this  problem. 

b.  Transfer  Function 

Let  the  input  to  the  system  be  a  plane  wave,  that  is.  g(t,r)  = 
_  J -xLiJ  Vo-')  wjtjt  freqUency  fo  and  spatial  frequency  vector  vQ.  The  output  of 
the  system  is  given  by  Eq.  (1)  as 


f(f-r) 


(-j‘2w(f0T-i/c,-T0) 


h(  r.rc,:r)  dr  d r0. 


(3a) 


<p(t  r) 


'jlnfj-vo-i) 


2*(f0r-v0 


•r0) 


h{r, r0;r)  dr  dr0, 


(3b) 


4 


or 


¥>(<,!>  =  h(r,r„;r)  } 


/=/,’ 

V=Vo 


(3c) 


where  $T  f represents  the  time-domain  Fourier  transform  and  ^  _v  represents  the 
spatial- domain  Fourier  transform.  The  subscripts  under  the  symbol  ^  represent  the 
variables  involved  in  each  transform.  We  define  H(fu\i),  the  system  transfer 
function,  as  follows: 


nOO  , 

r-j2ir(jr-wi c)  jT  ^ro  .  (4) 


Using  this  definition,  the  output  of  the  system  can  be  written  as 


(5) 


the  expected  result  from  linear  systems  theory. 

The  parameter  r0  in  the  integrations  in  Eqs.  (1)  and  (4)  represents  the 
vectorial  difference  between  the  observation  point  r  and  the  point  source  position  rs. 
The  integrations  must  be  performed  taking  r  as  a  constant,  that  is,  by  changing  the 
source  position  rs  so  that  r0  assumes  all  possible  values  in  Ft3,  and  the  results,  y?(t!,r) 
or  Hifo^oj),  are  valid  for  that  "fixed"  observation  point  r. 


c.  Response  to  a  Time-Harmonic  Point  Source 


When  evaluating  the  transfer  function,  it  is  convenient  to  use  a 
time-harmonic  point  source  as  the  input,  g(L,i)  =  6(r-rs)  with  frequency  jc 

and  location  rs.  The  output,  from  Eqs.  (1)  and  (4),  is  given  by 


<p(t,  r)  =  f  fc7"2^4  T)  6(r-rs-r0)  A(r,r0;r)  dr  dr0  , 
^  -  «*'  -  00 

p(t,r)  =  f  ft(r.r-rs;r)  dr. 


(6a) 


’Ur)  =  eJ'2T^/5r_r{  h ( r,r0;r)  } 


I  =  Jo 

rc=  r-rs 


(6b) 


(6c) 


or 


r(f-r)  =  Hifr.v -0:r)  )  .  (6d) 

d.  Tiansform  Relationships 

(1)  Fourier  Transform  Definitions  and  Notation.  The  time-domain 
Fourier  transform  of  a  function  g{  Lt)  =  g{  t.x.y.z )  and  its  inverse  ar^  given, 
respectively,  by 


G(l  r) 


$t_j{  ?Ur)  }  =  g(t.r) 


(7) 


9M  =  rt[j{  Gif, r)  }  =  J  G(/r) 

**  00 


(6) 


The  spatial  Fourier  transform  and  its  inverse  are  given,  respectively,  by 


G(L,v)  =  £r  |/  {  g{t,i)  } 


(9) 


and 


g{tr)  =  ^[v  {  G{t,v) }  =  j*  G{  t,v)  t~j^v'Tdv  , 


(10) 


where  dr  =  dx  dy  dz  ,  dv  =  df, c  d/.  d£  .  and  yr  =  xjk  +  y/y  +  cjz  .  The  above 
transforms  involve  triple  integrals  and  can  be  split  into  three  spatial  transforms  in  x 
(jk),  V  ( jy )  and  c  (j£).  For  example,  the  transforms  in  x  (j^)  are  written  as 


G{tS,y,z)  =  {  g(l,i,y,z)  }  =  J  g{t,x,y,z)  J “*kTdr  ,  (11) 


and 


g{t,x,y,:)  =  £  {  G{tS,y,z)  }  =  f  G{t,^,y,z)  e~*  “ThT<lf)L  .  (12) 


The  transforms  in  y  (fy)  and  2  (j£)  follow  the  same  notation.  The  full  space-time 
transforms  can  be  written  as 


G(M  =  {  9(0)  }  =  5,.^.^  Sz_k  {  9(0)  }  ,  (13a) 


or 


G{f,v)  =  f  f  g(tj)  dt  dr  , 

^  "  oc*'  “  00 

(13b) 

and 

S|,r)  =  <  ffu-)  1  =  S-.^  1  gi» } . 

(14a) 

or 

?Ur)  =  J  f  G(f.u)  dfdv. 

~  vC:  “  <X 

(14b) 

The  above  relationships  are  valid  fur  the  impulse  response/transfer  function 

is  replaced  by  t  and  r  is  replaced  by  r0  =  (it,Jfo>2o)>  that  is. 

pair  if  t 

Wf/rG;r)  =  ^r.y{  h(r;r0;r)  }  , 

(15) 

HyUr.r)  =  iNro_v{  h(r,r0;r)  }  , 

(16) 

(17) 


H{t,Lyo,Zo\r)  =  h(r,ib,yo,2o;r)  }  , 


and 

=  $T.fiXo.v{  MwO  }  .  (4) 

The  same  applies  for  the  respective  inverse  transforms. 

(2)  Input -Output  Relationships.  During  our  developments,  we  will 
need  to  express  the  output  of  the  system,  ^(t,r),  in  terms  of  the  transformed  input 
function  G(fv)  and  the  transfer  function.  From  the  Fourier  transform  definition 


g(t,x) 


ej'2r(ft-vr) 


G(fv)  df  dv  , 


(18) 


follows  the  transform  property 


g(t-r. r-r  ) 


J2x{f 


t-j2x(fr-vi0) 


G(fv)  df  dv  . 


(19) 


Substituting  this  expression  into  Eq.  (1)  yields 


-r)  e-j2r(jT-wi0) 


G(fv) 


*  A(r,r0;r)  df  dv  dr  dr0  , 


(20a) 
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p(f,r)  =  f  f  eJ  "ir^t~l/‘r^G(/,i/)f  f  A(r,r0;r)  e  v  r°^r  dr0  df  dv,  (20b) 

*  ~  OCT  ~  OC  ^  '  OC*'  * 

MW) 


or 


noo 

CO)  «0!^r)  J2'{Jt-‘"')dliv  ,  (20c) 


which  is  the  desired  relationship. 

Analogous  to  a  time -variant  system,  for  which  the  frequencies  at 
the  output  are  different  from  the  input,  a  space -variant  system  causes  the  output 
spatial  frequency  vector  to  be  different  from  the  input.  We  will  use  f)  to  represent 
the  output  spatial  frequency  vector.  T-,e  input -output  relationship  in  the 
transformed  form  can  be  obtained  by  subs  ituting  Eq.  (20c)  into  the  Fourier 
transform  definition 


{  p(t-r)  ) 


(21) 


Doing  so  yields 


oc  /•oc  /•  00  /•  IX 


*(M  = 


G(rj, v )  H(t),v,t)  v  dr)  dv  e  @  T^dtdr,  ( 


X  "X  'X  *  X 


4>(//?)  =  f  f  fc/ft/.i')  e  j  “^v  r(f  ei~T(f^td)>dT)drdv,  (22b) 

V  -  X*  -  oc*'  “X  '  "  X  / 
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*050  « 


G(r},i/)  H(%v\ r)  e'j2^I/‘^'r  6(/-»?)  dr]  )  rfr  di/,  (22c) 


*U0  «  J  f  GU,v)  drdv  .  (22d) 

-  (VI*'  -  fY\ 


2.  The  Acoustic  Channel 

As  long  as  we  deal  with  small -amplitude  acoustic  signals,  the  governing 
wave  equation  is  linear.  We  will  restrict  ourselves  to  time -invariant  problems, 
where  the  properties  of  the  medium  do  not  change  with  time,  and  the  source, 
receiver  and  the  boundaries  are  motionless. 

Figure  3  is  a  pictorial  representation  of  the  acoustic  channel.  A  source  or 
transmitter  rectangular  array  "couples"  an  electrical  signal  to  the  medium  and  a 
receiver  array  transforms  the  acoustic  signal  back  to  electrical  form.  Both  arrays 
are  assumed  to  act  as  linear  filters.  We  will  not  be  concerned  here  with  the 
generation  of  the  electrical  transmitted  signal,  assumed  given,  nor  with  the 
processing  done  on  the  received  electrical  signals  at  each  hydrophone  position.  The 
acoustic  channel  can  be  represented  by  the  block  diagram  shown  in  Fig.  1.  which 
can  be  broken  down  into  the  three  basic  components  sketched  in  Fig.  3:  transmitter 
array,  medium  and  receiver  array.  Figure  4  shows  these  components  in  block 
diagram  form.  [Ref.  5;pp.  3-4] 
a.  The  medium 

The  propagation  of  the  acoustic  signal  through  the  medium  is  governed 


by  the  wave  equation 


Fig.  4.  Block-Diagram  Representation  of  the  Acoustic  Channel. 


1 

V"y?U r)  -  ~tj — ^77  ^(t,r)  =  0(i,r)  (23) 

c^(r)  dt 

and  appropriate  boundary  conditions,  where  f(t,r)  is  the  velocity  potential 
(m~s  ),  g(t,i)  is  the  source  term,  which  represents  the  rate  of  mass  injection  into 
the  space  per  unit  mass  or  rate  at.  which  fluid  volume  is  added  per  unit  volume 
(  s  1)  [Ref.  5;p.  1].  These  two  terms  are  represented  in  Fig.  4  as  the  output  and 
input  to  the  acoustic  medium  block,  respectively.  The  solution  of  Eq.  (23)  is 
obtained  from  Eqs.  (1)  or  (20c)  as 


0 c 

X 


^M(r,r0;r)  dr  d ru  , 


(24) 


or 


fM  =  I  I  °U'V)  H*U'V'l)  ejMft~V'T)^dv  ■  (25) 
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The  transform  of  the  output  $(/$  is  a  function  of  the  spatial 
frequency  vector  /?,  which  is,  in  general,  different  from  v)  due  to  the  space -variant 
nature  of  the  medium.  The  physical  meaning,  in  terms  of  ray  acoustics,  is  that  a 
non -homogeneous  medium  (H  is  a  function  of  r)  causes  change  in  the  direction  of 
propagation  of  the  sound. 

The  input -output  relationship  in  the  transformed  form  can  be  obtained 
from  Eq.  (22d).  resulting  in 

noo 

G(f»)  H^:r)  -r  dr  du  (%) 


If  the  system  is  space -invariant,  an  unbounded  homogeneous  medium 
for  example,  then  H  is  independent  of  r  and  Eq.  (26)  reduces  to  the  expected 
expression  for  a  time-invariant,  space-invariant  linear  system: 


<*>(//* 


■f 


GW  H  if 


■l 


rj2r(i/-p)  -r 


dr  d  is 


27a) 


d’iffi) 


Gif v)  H^ifv)  Hv-fi)  d v  , 


(27b) 


*{fP)=G(fP)  H^fP). 


(27c) 
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b.  The  Transmitter 

The  transmitter  array,  the  first  block  in  Fig.  4,  is  characterized  by  a 
complex  aperture  function  Ax(/r)  (Amp’*  s~^)  relating  the  electrical  input  x(f,r) 
(Amp)  to  the  source  distribution  g(t,i)  (s'*)  through  the  expression  [Ref.  5:p.  30] 

G(/r)  =  X(fa)  AT{fj)  .  (28) 

O  1  I 

The  far -field  directivity  function  £>t(/i/)  (m  Amp’  s’1)  is  related  to  the  complex 
aperture  function  through  the  transform  [Ref.  5:p.  38] 

=  ^T.v{  -4T(/r)  }  .  (29) 

Using  Eq.  (29),  the  spatial  Fourier  transform  of  Eq.  (28)  is  given  by 


Oij.v)  =  X(  j.v)  *  Z>T(./>) . 


( 30a) 


or 


G(f.u) 


£ 


X(/a)  da , 


(30b) 


where  the  symbol  (*)  denotes  a  three-dimensional  convolution  with  respect  to  v. 
Equations  (30a)  and  (30b)  give  the  relationship  between  the  electrical  input  and  the 
source  of  the  acoustic  signal.  The  relationship  between  the  input  electrical  signal 
and  the  velocity  potential  (the  output  of  the  medium  block  of  Fig.  4)  is  obtained  by 
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substituting  Eq.  (30b)  into  Eqs.  (25)  or  (26).  Let  us  introduce  the  simplifying 
assumption  that  the  electrical  input  signal  is  independent  of  position,  that  is, 
x(t, r)  =  r(t)  and  X(f-v)  =  X(f)  $(i /)  .  This  is  justified  by  the  fact  that  the 
complex  weights  applied  to  the  electrical  signal  (array  shading  and  beam  steering) 
are  taken  into  account  in  Z?T(/i/),  the  fax  -field  directivity  function.  With  that 
assumption,  Eq.  (30a)  becomes 

G(U)  =  X(J)6(v)lDT(f,v),  (31a) 


or 


G{f.v)  =  X(j)  Dt(M 


(31b) 


which,  when  substituted  into  Eq.  (25),  gives  the  required  relationship  between  the 
input  electrical  signal  and  the  velocily  potential  at  the  output  of  the  medium  block, 


vHirl  =  J  j  A’I/IOjUiOH.U™) 

-  A* 


(32) 


A  lurther  restriction  we  impose  to  our  problem  is  that  the  transmitter  array  is 
composed  of  point  sources,  has  rectangular  shape,  is  parallel  to  the  XY  plane,  as 
shown  in  Fig.  3.  and  is  centered  at  position  rs  =  (j^,^,^).  The  array  has 
dimensions  Mt  x  Aft  (first  dimension  along  the  X  direction  and  second  dimension 
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along  Y),  with  interelement  spacing  dx t  meters  along  X  and  dyt  meters  along  the  Y 
direction.  Then,  the  far -field  directivity  function  for  At t  and  Nt  odd  is  given 
by  [Ref.  5;ch.  4] 


where  rpci  =  (pdxt  +  qdyt  +  ys,  2s)  represents  the  position  of  each  element  and 
rpq(/)  =  flpq (/)  e-7  is  the  complex  weight  associated  with  the  element  at 

position  rpq. 

c.  The  Receiver 

The  receiver  array,  the  last  block  in  Fig.  4.  is  characterized  bv  a 

complex  aperture  function  A  (/r)  (V  s  m’°)  relating  the  velocity  potential  tf(f.r) 

.0 

to  the  electrical  signal  y(t. r)  ( V  m  through  the  expression  [Ref.  5:p.  31] 


b'(/r)  =  $(/r)  AR(/r)  .  (34) 

As  stated  earlier,  we  are  not  concerned  with  the  signal  processing  done 

on  the  received  signal.  Our  purpose  is  to  predict  the  waveform  at  the  output  of 

each  receiver  array  element.  So,  we  consider  the  term  .4R(/r)  as  pertaining  to  a 

single  transducer  element.  The  quantity  y(U r)  is  the  electrical  signal  distribution 
-3 

(V  m  )  due  to  the  velocity  potential  <p(t,i)  at  a  point  r  on  the  transducer.  In  order 


IK 


to  obtain  the  total  electrical  signal  y(t)  at  the  transducer  terminals,  we  must 
integrate  y(t, r)  over  the  whole  transducer  volume,  that  is,  integrate  with  respect  to 

r, 

n/)=  f  nw*. 

'  oc 


n/i  = 


dr  ■ 


(35) 


by 


For  a  point  sensor  located  at  rr.  the  complex  aperture  function  is  given 


=  A(f)  6(T-zr)  ,  (36) 

which  results  in  an  electrical  signal  at  the  terminals  given  [  see  Eq.  (35)  ]  by 

Y(f)  =  <h(/rr)  A(f)  ■  (37) 

If  the  transducer  is  a  pressure  sensor,  then  .■!(/)  is  proportional  to  the  frequency  /, 
that  is,  A(f)~jf  and  Y(f)  ~  jfA>(J,Ts)  ~  pressure  .  For  simplicity,  we  take  A(f) 
as  unity,  and  a  pressure  sensor  can  be  simulated  by  inserting  a  suitable  filter  into 
the  block  diagram  of  Fig.  4.  So,  the  output  electrical  signal  will  be  numerically 
equal  to  the  velocity  potential.  The  receiver  array,  like  the  transmitter,  is  a 
rectangular  array  of  point  sensors  of  dimension  Mr  *  Nr  ,  parallel  to  the  XY  plane, 
with  interelement  spacings  dxr  and  dyr  meters,  centered  at  rr  =  (Jr,  J/r,  ^r)-  The 
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output  electrical  signal  at  each  element  ( m,n ),  given  by  Eq.  (37)  with  A(f)  =  1 
and  rr  =  rran  =  (to  <4r  +  ,  n  dyT  +  yT ,  Zr),  is 

W/) -*(£«).  (38) 


Vmn  (  ^  lp(  t,  rmn  )  • 


d.  Overall  Transfer  Function 

The  overall  transfer  function  relates  the  output  electrical  signal  ymn(t) 
at  each  point  sensor  (m,n)  of  the  receiver  array  to  the  electrical  input  signal  x(<)  at 
the  transmitter  array.  Substituting  Eq.  (39)  into  Eq.  (32)  we  obtain,  for  the  output 
signal, 


!/„(<)  =  f  °T a.  (40a) 

**  -  ft'1'  -  IV 


J00  C 00  'if 

W)J  CT(/»)  H,(/r;rmn) 


(40b) 


The  last  expression  can  be  rewritten  as 


,(<)-  J*.v(/)  Hao 


(40c) 
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or 


y^O)  -  rtlf{  XU)  HUi„) } ,  (40d) 

where  W(/rm)  is  the  overall  transfer  function,  given  by 

Joe 

Ot(/»)  .  (41) 

00 

Note  that  Eq.  (40c)  is  valid  only  when  the  input  electrical  signal  is  independent  of 
position  and  the  receiver  is  a  point  sensor  with  unit  frequency  response,  that  is, 
A(f)  =  1.  On  the  other  hand,  Eq.  (32)  assumes  only  the  independence  of  the  input 
electrical  signal  with  respect  to  position. 

B.  THE  RANGE  INDEPENDENT  A  COL  STIC  CHANNEL 

In  this  section  we  will  derive  the  expression  for  the  overall  transfer  function  for 
the  particular  case  of  a  iange  ir.,1  pendent  medium.  With  reference  to  Fig.  3,  the 
range  independent  medium  characteristics  change  only  with  the  Y  (depth) 
direction,  and  the  boundary  conditions  are  independent  of  j  and  z,  that  is,  the 
problem  presents  cylindrical  symmetry. 

1.  Medium  Transfer  Function 

In  order  to  obtain  the  medium  transfer  function,  we  will  solve  the  wave 
equation  using  a  time -harmonic  point  source,  that  is, 

1  £ 

V^Ur)  -  -77 - 7  y?(f,r)  =  <5(r-rs)  ,  (42) 

c-(  y)  dr 
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with  suitable  boundary  conditions  .  The  solution,  as  shown  in  Eq.  (6d),  is  of  the 
form  p(i, r)  =  £(r)  wjth  £(r)  satisfying  the  Helmholtz  wave  equation 

V2^(r)  +  k2(y)  £(r)  =  6(r-rs)  ,  (43) 

where  k(y)  =  ‘2rf  /c(y)  is  the  depth  dependent  wave  number.  Applying  the  spatial 
Fourier  transform  with  respect  to  xand  2  to  Eq.  (43),  we  obtain 

K(y)  H(jUO  +  ^H(J Uj»  =  «(»-*> .  (44) 

dy 

.0  .  .  0  .  .  .  O  .0. 

where  ky(y)  =  k~{y)-4x‘'(jx~+  jk“)  and  H(£,y,jk)  is  the  spatial  Fourier  transform 
of  £(r)  with  respect  to  x  and  2.  Solutions  to  Eq.  (44)  are  usually  obta.  d  by 
approximate  methods,  most  notably  the  WKB  approximation.  We  can  write  the 
solution  to  Eq.  (44),  without  loss  of  generality,  as 


=-U*,y,Jz)  =  o(  y) 


✓ 


(45) 


The  (transformed)  velocity  potential,  obtained  from  Eq.  (45)  and  the  assumed 
<p{t, r)  =  £(r)  is  given  by 

*UAv,0  =  H(A.!/.jt)^2’r/»'.  (46) 


oq 


The  velocity  potential  is  given  by  the  inverse  Fourier  transform  of  Eq.  (46) 

with  respect  to  land  z,  namely 


*>(1,  r)  = 


(47) 


Next,  we  will  relate  these  results  to  those  obtained  in  Section  II-A-lc,  the 
response  of  a  linear  time -invariant,  space -variant  system  to  a  time-harmonic  point 
source,  as  given  by  Eq.  (Gd ),  which  can  be  written  as 


A  nr. 


e'j  2*fyyoe-j  2  irfzZodl/ 


(48) 


where  r0  =  (^,  Jfc,  k>)  =  (j>£s,  y~ys,  z-*.)  and  dv  =  djx  dfy  d '£■  In  order  to 
compare  Eqs.  (47)  and  (4p),  we  rewrite  Eq.  (47)  using  the  term  defined  in 
Eq.  (45)  as  follows: 


'r  ( hr'  = 


*  JzZs)  p-j2irfxT  p-j'2tjzz 


dk  dk 


(4Qa) 


y{  t,  r)  = 


(49b) 


Comparing  the  Fourier  transform  expressions  in  Eqs.  (48)  and  (49b),  the  term 
obtained  from  the  solution  of  the  inhomogeneous  wave  equation  when  the  source  is  a 


unit  amplitude,  time-harmonic  point  source,  is  seen  to  be  related  to  the  medium 
transfer  function  through  the  transform 


=  ry'o.^{  ).  (50a) 

or 

niUo.kv)  =J  H,a k,h,kv)  ,  (50b) 

*  00 

where  y  has  been  substituted  for  r  as  an  argument  of  H  in  order  to  stress  the 
depth -only  dependence  of  the  medium  transfer  function.  Notice  also  that  the 
arguments  of  'l'  have  been  written  in  accordance  with  the  linear  systems  notation 
introduced  in  Section  II-A-1. 

The  transfer  function  can  then  be  obtained  as  the  direct  spatial  Fourier 
transform  of  with  respect  to  y0  =  y~ys.  In  the  next  section  we  will 

see  that  this  last  step  is  not  necessary,  so  long  as  the  transmitter  array  is  composed 
of  point  sources. 

2.  Output  Electrical  Signal 

The  overall  transfer  function  for  a  range  independent  medium  can  be 
computed  from  Eq.  (41)  by  substitution  of  the  expressions  for  and  the 

medium  transfer  function.  As  we  just  mentioned,  it  will  not  be  necessary  to 
compute  the  medium  transfer  function  if  the  transmitter  elements  are  point  sources. 
The  directivity  function  for  the  transmitter  array  is  a  sum  of  terms  of  the  form 
cpq(/)  r  ~  P9  ,  where  rpq  =  (xp  ,  yq  ,  2*)  is  the  position  of  each  transmitter 
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element.  Substituting  that  general  summation  term  into  Eq.  (41)  and  expanding 
the  exponential  factors  yields  the  result  that  the  overall  transfer  function  is  a 
summation — over  p  and  q — of  terms  of  the  form 


or 


'pq 


(/,f 


rJ 


r-J 24(#„-.v,)  t-j M<VV  HULMM  i» 


where  x^  =  xr+m  dxr  and  yn  =  yr+n  dyr  are  the  (x,y)  coordinates  of  the  receiver 
array  elements,  and  the  (x,p)  coordinates  of  the  transmitter  array  elements  are 
xp  =  rs+p  dxt  and  yq  =  ys+q  dyt  .  Recall  from  Section  II -A -1  that  the  parameter 
r0  is  the  vectorial  difference  between  observation  point  and  source  position.  The 
integral  in  this  last  expression  is  the  inverse  spatial  Fourier  transform  of  the  transfer 
function  Hu(  f,k-  fv-k'-V n)  with  respect  to  k  k  and  k  which  yields  a  function  of 

ll 

(Xj,  yc,.  Zr_. ),  where  x c  =  i^-xn  .  yt-,  =  yn -yn  ■  and  =  rr-Cs  .  But  the  transform 

r  ’-1 

with  respect  to  k  is  already  given  by  Eq.  (50b)  as  the  function  ^  (f,fx,yu-fi\y) 
obtained  from  the  solution  of  the  wave  equation.  As  a  result,  we  can  rewrite  the 
above  last  expression  as 


^if-Uy-ky) 


-J2*fz:o 


dk  dk  ■ 


Therefore,  the  overall  transfer  function  for  the  range -independent  medium  is  given 

t>y 


Hi  -1  hA 


A  00  J*  00 

^  X  cpq^J  I 

p  =  _m t-i  =_#t-i 


rj  2tjT; 


*  e  -/  — o  dfz  dfi  (51) 


The  output  electrical  signal  at  each  receiver  position  is  given  by  Eq.  (40c). 
repeated  here  for  convenience: 


wo-  f'w)  nutm) 


(40c) 


where  is  given  by  Eq.  (51).  Again,  we  stress  that  Eq.  (40c)  incorporates 

three  assumptions  regarding  the  linear  system: 

•  The  input  electrical  signal  is  independent  of  position  r,  that  is.  x(f,r)  =  x(  t)  . 

•  The  receiver  array  elements  are  point  sensors,  that  is.  A  (f, r)  =  ,4(/)  r-rr )  . 

•  The  sensors  have  unit  frequency  response,  that  is.  A(f)  =  1  . 

Equation  (51)  assumes,  in  addition: 

•  The  transmitter  array  elements  are  point  sources. 

•  The  medium  has  range  independent  characteristics  and  boundary  conditions. 

•  Transmitter  array  geometry  is  as  described  in  Section  I-A-2b. 


a.  An  Alternative  Representation 


We  will  now  compare  our  results,  as  given  by  Eqs.  (51)  and  (40c),  with 
the  results  obtained  by  solving  the  wave  equation  for  a  range -independent  medium 
in  cylindrical  coordinates  and  written  in  terms  of  a  Fourier-Bessel  integral. 
DiNapoli  and  Deavenport  [Ref.  4:p.  80],  for  example,  use  that  representation  as  a 
starting  point  for  the  description  of  numerical  methods  applied  to  range  independent 
problems. 

The  integral  in  Eq.  (51)  can  be  seen  as  the  overall  transfer  function  due 
to  a  single  point  source.  From  this  point  of  view,  is  just  the  weighted  sum 

of  the  transfer  functions  due  to  each  transmitter  array  element,  as  one  would  expect 
for  a  parallel  connection  of  linear  systems.  These  elementary  overall  transfer 
functions  are  given  by 


The  overall  transfer  function  ran  then  be  written  as 


Aft-1 
— T 


A-’t-l 

- 


M/rmn)  = 


r  _  M-i  -  _  Wt-i 
p  -  2“  9  ‘  — T 


£  W/)  Wpq(/rnJ  •  (W) 


Let  us  make  a  change  of  variables  in  the  integral  in  Eq.  (52),  exploiting 
the  fact  that  the  dependence  of  4*  (/£,  $/„,.£;  yn)  with  f*  and  k  occurs  due  to  the 
factor  ky(y)  =  ^"  ( i/) -4  -t  j^)  in  Eq.  (44).  that  is,  ^  is  a  function  of  (^“-r  k ^)- 
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Introduce  the  variables  fT  and  7  such  that 


k  =  /r  cos  7  and  fz  =  fT  sin  7  , 


(54) 


which  is  a  change  from  "rectangular"  coordinates  (kJz)  *°  *he  "p°^ar"  coordinates 
(jf.,7).  With  this  change  of  variables,  Eq.  (52)  becomes 


«pqar™>=[  *ufr,y°^)  t, ;|  (55) 

J  C\  J  ft 


where 


^{fJrVo'yn)  =  Vif-k-yrSWn) 


/x=/r  cos  Tr 

/z=/r  sin  7 


(56) 


Equation  (55)  can  be  written  as 


^pq(/rmn  ( 


J  T 

I  4M«  X  f  rJ-'/rflo  rf/r  , 

•'ll  •'o 


(5‘ 


where  sj/i  a  =  i0  /  R0  ,  cos  a  =  z0  j  R0  ,  and  /?0  =  v  Jq  +  zo  •  Because  of  the 
periodicity  of  the  yine  function,  the  term  o  can  be  dropped  from  the  exponent 
without  altering  the  result.  The  last  integral  can  then  be  recognized  as  the 
zero-order  Bessel  function  of  the  first  kind  [Ref.  6:p.  184],  that  is, 
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(58) 


u 


r2* 

z]  =  tA  6 

•'rt 


„ - jz  sin 


^  d') 


which  leads  us  to  the  Fourier -Bessel  integral  representation 


^(iu)  =  2 r  f  *(//r,y0;yn)  i0(2r/r/?0)  /r  dfT 

j  C\ 


(59) 


Equation  (59)  is  the  same  as  Eq.  (3.1)  in  [Ref.  4)  if  we  perform  the 
trivial  change  of  variables  kr  =  Jrrjj.  There,  'I'  is  interpreted  as  the  depth  dependent 
Green's  function.  Equation  (59)  is  the  starting  point  to  several  approaches  to  the 
range  independent  problem,  both  numerical  and  analytical  [Ref.  4:pp.  80-134]. 
Using  the  Fourier -Bessel  integral  representation,  the  overall  transfer  function 
H{  f. rnui )  can  be  written,  from  Eqs.  (53)  and  (59).  as 


A/t-1 

- T — 


A’t  -1 

- T- 


/*00 

n>  =  ^  '-pqi/)  -*1  ^(//r-J/oUVn)  /(,(27r/r/?0)  /r  d/r  .  (60) 


Mt-1  .Vt  - 1 

p  - - L__  y  = - 


C.  THE  LAYERED  WAVEGUIDE 

In  this  section  we  will  derive  the  function  for  the  case  of  a  three  fluid  medium 
waveguide,  commonly  referred  to  as  the  Pekeris  waveguide. 
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1.  Statement  of  The  Problem 

Figure  5  shows  the  relevant  physical  aspects  for  our  derivation.  The  point 
source  is  at  position  (rs,ys,2g)  in  medium  II  and  the  receiver  is  at  ( x,y,z )  in  the  same 
medium.  Medium  II,  with  characteristic  impedance  pjCj,  is  bounded  at  y  =  0  ,  the 
surface,  by  a  fluid  with  characteristic  impedance  plcl  and  at  y  =  D ,  the  bottom,  by 
another  fluid  with  characteristic  impedance  /?3c3.  As  described  in  Section  B-l,  in 
order  to  obtain  4*,  we  need  to  solve  the  Helmholtz  wave  equation,  Eq.  (44),  and 
isolate  the  terms  with  y  dependence  from  the  solution,  as  shown  in  Eq.  (45).  An 
inspection  of  those  two  expressions  show  us  that  by  setting  xs  and  zs  to  zero,  the 
solution  of  the  wave  equation  will  be  automatically  the  function  4'.  The  wave 
equations  that  describe  the  propagation  in  the  three  media  are: 

o 

kyj  4',  +  - ry  4>,  =  0  ,  (61) 

dy~ 

4>*2+^-2*2  =  <*(r3/s),  (62) 

dy 

and 

o 

*V3  ^3  +  ^3  =  ^  >  (63) 

ay 

o  o  o  o  o 

where  kyx  =  k\- 4a>'(jk“+  jk")  ,  kt  =  2*f /ci  and  the  subscripts  1,  2,  and  3  refer  to 
the  three  different  media.  The  argument  of  4*  has  been  dropped  for  simplicity.  We 
seek  the  solution  4>r  because  we  shall  be  concerned  only  with  those  problems  where 
both  source  and  receiver  are  in  medium  II.  The  boundary  conditions  are  the 
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Fig.  5.  The  Three  Layer  Waveguide. 
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=  A*  & '  y  +  R,  b  y.y,<y<D, 


(68) 


and 


#3  =  Bz  e~*y*  y  ,y>D.  (69) 

The  solution  in  medium  II  has  been  split  in  two  in  order  to  account  lor  the  point 
source,  as  usual  in  the  Green's  Function  determination.  Two  more  boundary 
conditions  are  necessary: 

•  Continuity  of  the  solution  at  V  —  V*- 

•  Discontinuity  of  the  first  derivative  of  at  y=  j/s. 

The  value  of  the  discontinuity  can  be  found  by  integrating  Eq.  (62)  over  a  region 
|  y-ys\  <  f  and  taking  the  limit  as  e  -*  0  .  Performing  this  operation  and  taking  into 
account  that,  in  the  limit,  the  integration  of  the  continuous  function  is  zero,  we 
obtain 


\k2  dy  =  1  , 


(70a) 


or 


(70b) 
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a.  Boundary  Condition  at  y  =  0 


In  order  to  obtain  the  expression  that  represents  the  continuity  of 
pressure  at  y  =  0,  we  equate  the  value  of  pressures  due  to  ty,  and  vI'2a.  Substituting 
the  assumed  solutions — Eqs.  (66)  and  (67) — into  Eq.  (64)  yields 


Pi  —  Pi  (c42a  +  B2a)  . 


(71) 


For  the  normal  component  of  particle  velocity,  Eq.  (65)  and  the  assumed  solutions 
yield 


K\  -41  =  *v2  (Aa-B2a>  • 


(72) 


Substituting  Eq.  (71)  into  (72).  is  eliminated  and  the  relationship  between  the 
amplitudes  of  t,he  solution  4/oa  is  found  to  be 


—  =  /£> 

A,. 


2  k 
po  y*- 


-  k 


V  l 


plk 

'•..O 

1‘r,  >• 


(73) 


When  k„ 2  is  real,  corresponding  to  propagating  waves  in  medium  II,  the  factor  ft,,  is 
physically  interpreted  as  the  plane  wave  reflection  coefficient  at  the  surface. 

b.  Boundary  Condition  at  y  =  D 

Following  the  steps  of  the  previous  section,  using  the  forms  of  the 
assumed  solutions  for  4*3  and  at  y  =  D  .  we  find 
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P's  -4S  ^jk^D  =  Pi  {A2b  A* D  +  e AsD  )  ,  (74) 

and 

-*y3  ^3  eA3D  =  *y2  (A,b  A2D-B2b  fAD  )  .  (75) 

Substituting  Eq.  (74)  into  (75),  .43  is  eliminated  and  the  relationship  between  the 
other  two  amplitude  factors  is  found  to  be 


•4ob 

— -  e 

B2h  ‘ 


- 1 

pi 


ft . 

—  k 

Pi  y 2 


y3 

—  e-j'2kyiD 


\V3 


(76) 


When  ky,,  is  real.  A;,  has  the  physical  meaning  of  the  plane  wave  reflection 
coefficient  at  the  bottom. 

c.  Source  Condition 

Using  the  results  of  the  two  last  sections.  Eqs.  (73)  and  (7b).  the 
solutions  in  medium  II  become 

=  42a  ( A’?y  +  A,  c  ),  0  <  y  <  ys  .  (77) 

and 

*2b  =  B?h  (Hi*  t~J  2kY?DAiy  *  c  A^  ),  ys  <  y  <  D  .  (78) 
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The  condition  of  continuity  of  at  y  =  ys  ,  applied  to  the  above  two  equations, 
gives  the  relation  between  the  amplitudes  ^42a  and 


A,a  «23  e~j2ky*D  eA2  y*  +  e-tty*  h 
B2b  As  y»  +  Ril  e"A2  Vs 


(79) 


Finally,  apply  the  condition  on  the  derivative  of  ^  at  y  =  ys  ,  as  given  by 
Eq.  (70b),  to  obtain: 

B2b  (R23  e~J  2kytD  A2^-e  A2^  )-^2a  (  A2^,  f^y2!/s  )  =  —  ,  (80) 

jky2 


which,  together  with  Eq.  (79),  forms  the  necessary  system  of  equations  to  obtain  the 
two  amplitude  factors.  B2b  and  .42a.  Upon  substitution  of  the  expressions  for  R>b 
and  ,42a  into  Eqs.  (77)  and  (78),  we  obtain  the  solutions  ^2a  and  ,  which  form 
the  transformed  transfer  function  'M/JU J/&,i;?/)  we  seek.  Written  in  an 
abbreviated  form,  the  final  expression  is 


i  H23  f'i2iUD  A2 

*Ujfc,*,jfc:y)  =  A - 


2Jk 


V2 


,  p  n  -j  2k 
1  -  fu  i  R23  t  > - 


(As  y<  +  Ha,  eA2  K  )  , 


(81) 


where  0<t/<D,0<ys<D,  and,  as  usual,  t/>  is  the  greater  and  y<  the  smaller 
between  y  and  ys  . 


3.  Special  Cases 


In  order  to  check  the  solution  for  the  layered  waveguide,  we  will  now 
derive  the  expression  for  the  output  electrical  signal  in  two  special  cases  whose 
analytical  solutions  are  known.  These  results  will  be  used  to  check  the  computer 
implementation.  In  both  cases,  the  input  will  be  a  time-harmonic  point  source. 
cl  Unbounded  Homogeneous  Medium  (Free  Space) 

An  unbounded  homogeneous  medium  can  be  represented  by  Fig.  5  if 
the  sound  speeds  c{  and  densities  are  equal  for  the  three  media.  In  this  case,  the 
reflection  coefficients,  as  given  by  Eqs.  (73)  and  (76).  become  zero.  The 
transformed  transfer  function — Eq.  (81) — reduces  to 

=  j—  e~jkyty>  e>ky*y<  .  (82a) 

-ky  O 

=  j—  r^y>  -  y<) .  (82b) 

2  *,.<> 


or 


^(//x.  &  ■£••?/)  =  j— -  C~JkY'l\y  =  j— - 


2  A 


v- 


">A 

-Ay  2 


(82c) 


Note  that  the  right  hand  side  of  Eq.  (S2c)  is  expressed  in  terms  of  y0  only,  meaning 
that  the  unbounded  homogeneous  medium  is  space -invariant,  as  expected. 

The  overall  transfer  function  H(Ji)  for  the  case  of  a  point  source  is 
given  by  Eq.  (51)  with  Mt  =  Nt  =  c, ,(/)=  1.  Substituting  Eq.  (82c)  into  Eq.  (51) 
yields 


The  above  integral  is  the  expansion  of  a  spherical  wave  into  plane  waves  [Ref.  7:p. 
B-5]  and  the  resultant  overall  transfer  function  is 

-jknR 

W(/r)  =  -  * - ,  (84) 

4tR 

r-n - 7 - 7  / - n - 7 - 7 

where  R  =  ||  r-rs  ||  =  y/  +  y“  -  =  >/( x-xt)“+{y-yt)~+{z-it)~- 

The  output  electrical  signal  y(t) — or  velocity  potential  ip(t, r),  in  our 
case — is  given  by  Eq.  (40c).  Substituting  Eq.  (84)  and  X(/)  =  S(f-fj) — for  the 
time-harmonic  source — into  Eq.  (40c)  yields 


1 

y(t)  - - 

4  IT  R 


J2rfj  -  k0R) 

C  > 


(65) 


where  h  =  2 x/.,  j c,,  .  Equation  (85)  is  just  the  expression  for  a  sphericaJ  wave,  since 
R  is  the  distance  between  source  and  receiver. 

b.  Surface  Reflection 

Consider  the  case  where  media  II  and  III  are  the  same,  that  is.  c2  =  c3 
and  Pi  =  p'i  ■  As  a  result.  R^  is  zero  and  ^  reduces  to 

*(f,U knky)  = y>  {^y<  +  «2,  f'^y2jK  ) :  (86a) 

2Ay2 
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*U.U %,k,y)  =  J—  (e'My>  - y<)  +  fl,,  e'-M15,  *y<)),  (86b) 

2kyl 

or 

*U,Ly°,ky)  =  J  —  ■  *1  +  ft,,  e-^ly  +  *‘l  ) ,  0  <  y,  0  <  y,  .  (86c) 

2iy2 

Comparing  Eqs.  (62c)  and  (86c),  we  see  that  this  problem  can  be  treated  by 
considering  a  homogeneous  medium  with  two  point  sources,  one  located  at  ys  with 
unit  strength,  and  the  second  located  at  -ys  with  strength  FU..  The  output 
electrical  signal  due  to  the  first  term  in  Eq.  (86c)  has  already  been  computed  in  the 
previous  section,  and  is  given  by  Eq.  (85).  In  the  second  term,  we  have  the  factor 
fio,.  also  a  function  of  fx  and  /z  through  the  dependence  on  kyi  and  k  When 
performing  the  inverse  spatial  Fourier  transform  to  obtain  F/(/.r),  it  would  be  easier 
if  /?21  were  constant.  If  we  consider  the  boundary  I -II  as  pressure  release  ,  that  is. 
pjp*  =  0,  then,  by  Eq.  (73),  FU{  =  -1.  In  this  case,  Eq.  (86c)  reduces  to 

♦(jU-jOw)  -  J—  .  ,87) 

O  l 

-  A  yO 

where  y £  =  y-(  ~ys)  ■  From  the  results  in  the  previous  section,  the  output  electrical 
signal  is,  by  superposition, 

y(t)  =  -  —  +  'ftlJot-KR')  (  (88) 

4j  R  4  icR' 
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where  R '  =  ||  r-r^  ||  =  Jt2Q  +  y'2  +  zj  =  y/{x-xs)2+{i ny^^z-z^)2  .  If  the  distance 
R  is  large  compared  to  the  depth  of  the  source  ys ,  a  simplification  is  possible.  From 
Fig.  6,  for  large  R  ,  Rx  RS-AR  ,  R'  x  Rs  +  AR  and  A R  x  y&  sin6  =  ys  yj  Rs  , 
where  Rs  is  the  distance  from  the  receiver  to  the  midpoint  between  the  source  and 
its  image.  With  these  approximations  and  letting  R  x  R'  x  Rs  in  the  amplitude 
factor,  Eq.  (88)  becomes  [Ref.  7:pp.  170,  408] 


y{t)  = 


j—  JWo*  -  koRs)  si 

2*RS 


ko  yS  Vr 


Rc 


(89) 


D.  SUMMARY 

The  transformed  transfer  function  $  {Ux.VoSw)  for  a  range -independent 
medium  is  the  solution  to  the  inhomogeneous  Helmholtz  wave  equation 

o 

K{y)  ^  #  =  Hy-ys) ,  (90) 

dy~ 

o  o  9  o  o 

with  suitable  boundary  conditions,  where  ky  =  A7-4t  (,&%  &*)  and  k  =  2* f /c. 

Under  the  conditions  stated  in  Section  II -B -2,  related  to  Eq.  (51),  the  overall 
transfer  function  for  a  range-independent,  medium  is  given  by 

Mt-l  Nt-1 

— T  ““ 2" 

««»-,)=  £  £  cpq(/)[  f  *(/£.*.&*.) 

Mt-l  Nt-1 

P  =  --V  q  =  -~V 

•  r’WAdjtii,  (5i) 
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where  the  transmitter  array  is  rectangular,  centered  at  (xs,ys,2s),  parallel  to  the  XY 
plane,  with  odd  dimensions  Mt  *  Nt  ,  complex  weights  c  (/),  and  with 

r  t 

elements  positioned  at  xp  =  is+px  ,  yq  =  y6+qy  ,  and  2  =  2S  .  Also,  x0  =  rm-rp  , 
y0  =  yn-yq  ,  and  z0  =  zT-zs  .  The  overall  transfer  function  /■/(/ r^)  can  also  be 
expressed  as 


Aft -1 


N,  -1 


r _ =  Jf 


X  X  CP q^f  ^UJv^y n) 

-  .Mid  -  .  M-l 


A/t-1  .  _  -  . 

~T  q  ~  7 


*  J0(2rJrRo)  fT  dfr ,  (60) 


j — - 57  /  7 - 7 

where  /r  =  /  -  /2“  .  /?,,  =  7  r"  -  2“  ,  and 


^(//r^o^n)  =  V  (f,hyo.fi-yn) 


/x=/r  7 

/z=/r  stn  7 


(56) 


The  output  elect riciil  signal  ymn ( /),  under  the  same  conditions  stated  for  the 
overall  transfer  function,  is  given  by 


,V)  =  j'’xu)HU*m,)'i'2wft<‘f, 


(40c) 


where  /f(/rran)  is  given  by  Eqs.  (51)  or  (60). 


4'J 


The  layered  waveguide  illustrated  in  Fig.  5  is  characterized  by  the  transformed 
medium  transfer  function 


.  a.  e-i2k)2De}^2  v>  + 

b,A;y)  =  j— - : - 

2ky 2  l-  tf2I  tf23  e~i  2ky*D 

,  (A>  X  +  R,, '-***<)  ,  (81) 


which  is  valid  for  0  <  y  <  D  ,  0  <  ys  <  D .  and  where  fi2|  *s  given  by  Eq.  (73)  and  /ko 
by  Eq.  (7b).  If  all  three  media  have  the  same  acoustic  characteristics,  then  Eq.  (81) 
reduces  to 


^U-jx-yu-Sw) 


and  the  corresponding  output  electrical  signal  is  given  by 


1 

y(t)  = - 

4  rB 


£>-*£*  -  k,,R) 


(82c) 


(85) 


where  R  is  the  distance  from  source  to  receiver.  In  the  case  where  media  II  and  III 
are  the  same  and  the  surface  is  pressure  release,  the  function  is  given  by 


y(lh  yo-k-y)  =  j-^—  (f  jk>M-e  ikydy°\) ,  (87) 

2 


43 


where  y'0  =  y  +  y6  .  The  output  electrical  signal,  when  the  source -receiver  distance 
is  large  compared  to  the  source  depth,  is  given  by 


y(  0  = 


.j—  jwj-w 

2*Rt 


ko  ys  yr 


R 


S 


(89) 


where  Rs  is  the  distance  from  the  receiver  to  the  midpoint  between  the  source  and 
its  surface  image. 
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III.  IMPLEMENTATION  AND  RESULTS 


In  this  chapter,  a  computer  implementation  of  the  range  independent  medium 
equations  is  described.  Although  only  the  medium  transfer  function  of  the  layered 
waveguide  is  implemented,  the  modular  nature  of  the  linear  systems  approach 
enables  us  to  analyze  the  results  and  draw  conclusions  valid  for  the  more  general 
case.  Equation  (60)  illustrates  this  point.  The  summations  and  the  coefficients 
c  (/)  in  that  equation  represent  the  transmitter  array.  The  Bessel  function 
J0(2~jrR0)  is  the  range  dependent  factor.  The  function  (f,fx,y0,fz]y) — from  now 
on  referred  to  as  the  medium  transfer  function,  for  simplicity — is  the  (range 
independent)  medium  dependent  factor. 

Our  purpose  is  to  have  a  working  algorithm  to  test  the  results  from  the  previous 
chapter.  Care  has  been  taken  to  enable  the  implementation  of  other  medium 
transfer  functions,  allowing  for  depth  dependent  speed  of  sound,  so  that  the  program 
can  be  used  as  an  academic  tool.  Speed  and  style  were  secondary  objectives. 

The  implementation  is  based  on  Eq.  (60),  the  Fourier -Bessel  integral 
representation. 

A.  IMPLEMENTED  EQUATIONS 
1.  Medium  Transfer  Function 

The  medium  transfer  function  $ (fjr,y0\y)  is  an  oscillatory  or  exponential 

2  2 

function  of  kyZ — as  seen  in  Eq.  (81) — which,  in  turn,  is  a  function  of  (fx  +  Jf  )  = 
o 

=  fr~.  The  function  ^  is  easier  analyzed  in  terms  of  ky7  than  in  terms  of  /r  = 
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=  yj  (/ / c^f  -  {kyi/2r)“  ,  which  is  nonlinear.  Therefore,  will  be  treated  as  a 
function  of  ky2,  and  the  integral  term  in  the  expression  for  the  overall  transfer 
function,  Eq.  (60).  will  be  modified  accordingly. 

a.  as  a  Function  ofky 2 

The  function  $  is  already  expressed  in  terms  of  fcy2  in  Eq.  (81).  The 
reflection  coefficients,  on  the  other  hand,  are  expressed  in  terms  of  lfcyl,  Jfcy2,  and  iy3. 
In  this  section  we  establish  the  relationships  between  ky2  and  the  quantities  jv  Jk  ,, 
and  kyz  .  These  relationships  will  enable  us  to  represent  the  reflection  coefficients  in 
terms  of  the  variable  A:y2,  as  well  as  to  change  the  variable  of  integration  in  the 
expression  for  the  overall  transfer  function  to  ky2. 

The  correspondence  between  ky2  and  /r  is  one-to-one.  When  /r  <  (//c2)  , 

ky2  is  taken  as  the  positive  root  (//e^-  ,  a  consequence  of  our  choice  of 

solutions,  Eqs.  (67)  and  (68).  When  fT  >  (// Cj),  ky2  is  imaginary,  corresponding  to 
an  exponentially  decaying  $.  When  all  factors  in  Eq.  (81)  are  multiplied  out,  is 
seen  to  be  composed  of  a  sum  of  factors  of  the  form 

A  e~fiy?a  , 

where  a  is  a  nonnegative  real  number.  In  order  to  have  t  be  a  decreasing 
exponential,  the  negative  root  must  be  chosen  when  ky 2  is  imaginary.  Similar 
reasoning  yields  the  same  results  regarding  Ay,  and  k  3,  which  cam  be  summarized  as 


follows: 


*yi  = 


-  *5  ,  *r  <  *i 


(91) 


-J 


*?  -  *?  .  *,  >  *i 


where  i  =  1,  2  or  3  ,  &r  =  2rfT  and  kx  =  2rj  fcy.  Using  Eq.  (91),  the  wavenumber 
y-components  kyi  and  iy8  can  be  expressed  in  terms  of  Jfcy2  as 


kyp  ~ 


£  -  £  +  *y2  » 


(92) 


2  2  2 
^2  2  i  >  K 


I — ^ - 7“ 

where  p  =  1  or  3,  and  k{  =  ^  —  Ay2  .  Using  Eq.  (91),  the  integrations  over  f( 

can  be  transformed  into  integrations  over  Av2.  Using  Eq.  (92),  the  reflection 
coefficients  /22,  and  /?28  as  given  by  Eqs.  (73)  and  (76),  can  be  expressed  in  terms 
of  *y3. 

b.  Air,  Water,  Fast  Bottom  Waveguide 

The  present  implementation  of  the  medium  transfer  function  assumes 
that  c(<  Co.  When  the  surface  is  almost  pressure  release,  which  is  the  case  of  the 
air-water  interface,  and  the  bottom  is  fast,  c3>  Cj  ,  the  medium  transfer  function  is 
composed  of  many  sharp  peaks,  which  presents  some  difficulties  for  integration 
routines. 

(1)  Poles  of  the  Medium  Transfer  Function.  When  the  surface  is 
pressure  release,  that  is,  /?2|  =  -1  ,  Eq.  (cl)  reduces  to 
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,  A>3  t^2ky^D  Jky^>+ 

^U^iVo\y)  =  -  - - - - sin(  k  2  y<  )  ,  (93) 

ky*  1+  #2  3  e‘J  2ky2D 


which  is  valid  for  0  <  y  <  D  and  0  <  ys  <  D  .  If  the  bottom  is  fast,  £,  >  *g  ,  and, 
from  Eq.  (91),  there  is  a  range  of  values  of  kr  for  which  Ay2  is  real  while  jfcy3  is 

imaginary ,  that  is,  k3  <  kr  <  k%  or  0  <  iy2  <  J  i2  -  k^  .  For  this  range  of  values 
of  ky2,  the  reflection  coefficient  H23)  given  by 


Pz 

ft  V  —  *y3 

ft3=  ^ - ,  (94) 

—  *y2  +  *y3 


becomes  a  unit  magnitude  complex  number, 


(95) 


where 


Pz  ] m{  ky3  } 
tan  i>  — - : - 

Pz 


Pz 

Pz 


(96) 


In  this  case,  the  denominator  of  the  expression  for  'I',  Eq.  (93),  may  have  zeros, 
given  by 


l  *  si,.-; 2i*i>_o, 


(97a) 
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i  +  e-i^D  -t)  =  o. 


(97b) 


Therefore, 


4>  -  kv2D  =  (2n  +  1)-,  n*  0,1,  2,...  (98a) 

y  2 


or 


where  0  —  ky2D  and  a  =  D  \  k2  -  k^  .  Solutions  to  the  above  transcendental 
equation  in  terms  of  ky2  =  0/ D  are  the  poles  of  the  medium  transfer  function,  for 
the  case  of  a  pressure  release  surface.  Eq.  (93),  and  a  fast  bottom.  The  physical 
interpretation  of  these  poles  is  that  they  represent  the  discrete  eigenvalues  ky2n 
corresponding  to  the  trapped  or  normal  modes  in  the  waveguide[Ref.  8:pp.  430-440]. 
The  graphical  solution  of  Eq.  (99)  [Ref.  S'.p.  438]  reveals  that: 
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•  Poles  exist  only  if  a  >  */2.  If  a>  t  ,  the  poles  will  be  separated  by  AJ ky2»  *(D. 

•  The  total  number  of  poles  is  given  by  the  integer  number  equal  to  l(o/x)-  0.5] 
or,  if  the  result  is  not  integer,  its  nearest  larger  integer. 


When  the  surface  is  almost  pressure  release,  that  is,  R2l  s  -1  , 
the  medium  transfer  function  presents  peaks  at  values  of  &y2  corresponding  to  the 
poles  computed  as  above. 

2.  Overall  Transfer  Function 

In  order  to  change  the  variable  of  integration  in  Eq.  (60),  the  interval  of 
integration  /r  €  (  0  ,  oo  )  is  subdivided  into  (  0  ,  / / )  and  (  / /c^  oo  ).  After  the 
change  of  variables  according  to  Eq.  (91)  and  a  rearrangement,  Eq.  (60)  becomes 


mn)  ~  I  2* Cp ql/J  ^  (MU  *y2  <^y2 

*T 


where  the  path  of  integration  P  consists  of  the  imaginary  axis  over  the  interval 
( -j  oc  ,  y  0)  plus  the  interval  (0  ,  f  / c)  on  the  real  axis.  The  summations  over  p  and 
q  are  as  indicated  in  Eq.  (60). 

B.  THE  SUBROUTINES 

In  this  section  we  describe  the  subroutines  and  present  the  results  for  some 
examples.  Also,  a  description  of  the  main  program  and  some  auxiliary  subroutines 
is  given,  which  will  help  to  clarify  some  aspects  of  the  implementation'. 


'The  main  program  and  auxiliary  routines  were  written  by  Dr.  Lawrence  J. 
Ziomek,  at  the  Naval  Postgraduate  School,  as  part  of  an  ongoing  research  project  in 
model -based  signal  processing. 


1.  Main  Program  and  Auxiliary  Routines 

The  main  program's  simplified  pseudo-code  is  shown  in  Fig.  7.  It  consists 
basically  of  calls  to  subroutines.  The  readin  subroutine  reads  the  input  data, 
consisting  of  transmitted  signal  characteristics,  transmitter  and  receiver  array  data, 
medium  characteristics,  and  control  data,  including  the  ocean  medium  transfer 
function  (OCNTF)  variable  used  to  select  a  transfer  function,  when  more  than  one 
type  of  transfer  function  is  implemented.  The  ccq  subroutine  generates  the  input 
electrical  signal  x(t).  The  chfmn  subroutine  computes  the  overall  transfer  function, 
Eq.  (100).  and  its  estimated  error.  The  phjrnn  subroutine  provides  a  tabular  output 
for  the  values  of  the  overall  transfer  function  and  its  error.  The  calycc  subroutine 
computes  the  output  electrical  signal,  according  to  Eq.  (40c). 


Main  Program  /computes  output  waveform  of  an  acoustic  channel/ 
EXECUTE  Readin  /read  input  data/ 

EXECUTE  Ccq  /generate  input  electrical  signal/ 

EXECUTE  Chfmn  /compute  overall  transfer  function/ 

IF  (flag  print  is  TRUE)  THEN 

EXECUTE  Phfmn  /print  overall  transfer  function/ 

END  IF 

EXECUTE  Calyce  /compute  output  electrical  signal/ 

End  Main  Program 

Fig.  7.  Simplified  Pseudo-Code  of  the  Main  Program. 


a.  Subroutine  Ccq 

The  subroutine  ccq  generates  the  input  electrical  signal  i(i),  using  the 
complex  envelope  representation.  From  the  specifications  read  as  input  data,  this 
subroutine  computes  the  time  samples  of  a  reference  complex  envelope  and  its 
Fourier  series  coefficients.  The  Fourier  series  is  truncated,  yielding  the  complex 
envelope  of  the  transmitted  signed. 

(1)  Complex  Signal  Representation.  An  amplitude  and  angle 
modulated  carrier  r(f)  can  be  represented  by  its  complex  envelope  f(t)  as  [Ref. 

5:pp.  177.  182] 


r(t)  =  Re{  f(t)  } 


(101) 


where 


U0  =  <■(!) 


(102) 


a(t )  is  the  (real)  amplitude  modulating  signal  (Amp),  9(t)  is  the  angle  modulating 
signal  (rad),  and  jc  is  the  carrier  frequency  (Hz).  The  Fourier  transform  of  the 
signal  is  related  to  its  complex  envelope  transform  by  [Ref  5:p.  187] 


R(f)  =  0.5[/?(/-  fc)  +  /?*(-/-  fz)]  . 


(103) 


For  a  rectangular -envelope  LFM  pulse,  a(t)  and  #(t)  are  given  by 


a(t)  = 


A  ,|«|  <  Tp/2 
0  , \t\>  Tp/2 


(104) 
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and 


0(t)  =  ^!—t2  =  Dpt2,  (105) 

T? 

where  A  is  the  (constant)  amplitude  (Amp),  T?  is  the  pulse  length  (s),  SW  is  the 

-2 

swept  bandwidth  (Hz)  and  Dp  is  the  phase-deviation  constant  (rad  s  ).  Both  S w 
and  Dp  are  allowed  to  be  either  positive  or  negative,  corresponding  to  the  up  chirp 
and  down  chirp  LFM  pulses,  respectively.  If  the  modulated  signal  r(t)  is  made 
periodic  with  period  T0,  then  the  complex  envelope  can  be  represented  by  the 
Fourier  series 


r(t)  = 


/ 


2 

k  =  -/ 


(106) 


where 


b 


k 


r(t)  e 


dt . 


1 1 07) 


and  J0  =  1/T0  .  Equation  (106)  assumes  a  complex  envelope  bandlimited  to  h’/Q  . 
From  Equation  (106),  the  transform  of  the  complex  envelope  is  given  by 


i 

R(f)=  £  (108) 

k  -  -f 
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If  the  complex  envelope  is  sampled  with  sampling  period  Ts  =  T0/(2K  +  1),  which 
is  basically  the  Nyquist  rate,  then,  from  Eqs.  (106)  and  (107), 


r(lTs)=  j  \J2*kllL, 
k=  -/ 


(109) 


L  )  r 


u  l  = 


(110) 


where  Ts  =  T0/L  is  the  sampling  period,  and  L  —  2 AVI  is  the  total  (odd) 
number  of  time  and  frequency  samples.  From  Eqs.  (103)  and  (108),  the  signal 
transform  is  given  by 


R(f)  =  0.5  £  {  K  Hf-  k  -  */<.) 
;•  =  -/ 


-  (in) 


from  which  the  sampled  signal  can  be  derived  as 


r(/Ts)  =  Re  |  £  bk  j2tkllL 

l  k  =  -/ 


(112) 


(2)  Subroutine  Description.  Figure  8  lists  the  pseudo-code  of  the 
subroutine  ccq.  The  main  input  data  for  the  subroutine  are: 


Pulse  length  Tp,  and  pulse  period  (or  pulse  repetition  period)  T0. 

Swept  bandwidth  Sk  and  pulse  amplitude  A. 

Maximum  order  of  frequency  component  to  be  processed,  KaiX.  The  number  of 
frequency  components  to  be  processed  is  (2  Kaix  +  1). 

Sampling  rate  multiplier  5k,  a  term  to  be  multiplied  by  the  Nyquist  rate  in 
order  to  obtain  the  actual  sampling  rate. 

Carrier  frequency  jc  . 

A  flag  (es/used  to  indicate  option  for  plotting. 


Subroutine  Ccq 

Compute  data  for  reference  signal: 

total  number  of  samples  L  for  the  complex  envelope  /Eq.  (114)/ 
sample  period  IS  for  the  complex  envelope  /Eq.  (113)/ 
phase— deviation  constant  DP  /Eq.  (105)/ 

Generate  samples  of  ref.  signal  complex  envelope  /Eq.  (102)/ 
Compute  complex  Fourier  coefficients  /Eq.  (110)/ 

IF  (flag  t'*t  is  TEUE)  THEN 

Compute  and  plot  transmitted  signal  /Eq.  (116)/ 

END  IF 

End  Subroutine  Ccq 


Fig.  8.  Simplified  Pseudo-Code  for  the  Signal  Generator  Subroutine. 


The  total  number  of  samples  for  the  complex  envelope  is  given  by 


L  -  To/Ts ,  (113) 

where  T0  and  Ts  are  the  pulse  repetition  period  and  complex  envelope  sampling 
period,  respectively.  The  sampling  frequency  fs  =  1/TS  is  the  product  of  the 
sampling  rate  multiplier  Sr  by  the  Nyquist  rate,  twice  the  complex  envelope 
bandwidth.  In  terms  of  the  input  parameters  listed  earlier,  Eq.  (113)  can  be 
rewritten  as 


h 

- - S 

LzTqSr  2  (  |5w  |  +  2/Tp  )  ,  (114) 

' - - - ' 

z  envelope) 

which  is  the  expression  used  in  the  subroutine.  If  the  number  of  samples  is  less  than 
the  number  of  frequency  components  to  be  processed,  (2  A'n,ax  +  1),  L  is  increased 
ordingly.  If  the  computed  value  of  L  is  even,  the  result  is  increased  by  one.  The 
phase-deviation  constant  Dr  and  the  sampling  period  Ts  are  computed  according  to 
Eqs.  (105)  and  (113),  respectively. 

The  subroutine  proceeds  to  compute  the  samples  of  the  complex 
envelope  of  the  reference  signal — using  Eqs.  (102),  (104)  and  (105) — at  the  time 
instants  t  =  l  Ts  .  The  Fourier  coefficients  are  computed  next,  using  Eq.  (110). 

The  signal  to  be  transmitted  is  defined  by  multiplying  the  Fourier 
coefficient  series  {6k} — computed  for  the  reference  signal — by  a  window  of  length 
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(2  A'max  +  1)-  The  simple  truncation,  which  corresponds  to  a  rectangular  window,  is 
avoided  because  of  the  ringing  it  would  cause  on  the  waveform.  Instead,  the 
Hamming  window 


w(n)  =  0.54  +  0.46  cos( t  n/KmiX)  ,  ®  -  ^max  >  (115) 


is  used.  Therefore,  the  transmitted  signal  and  its  complex  envelope  are  given  by 


x(/rs)  =  R‘ 


-K<  /<  A', 


(116) 


and 


fmax 

x(/rs)=  y  tflfjfc)  bkfP2*kllL  ,  -A'</<  A'  ,  (117) 

^  =  ax 

The  transform  of  the  transmitted  signal  is  given,  from  Eq.  (116),  by 

fmax 

A'(/)  =  0.s£  «ik){bkSU-Jc-kl)*bmt6(-}-l-k/J).  (lie,) 

^  =  ~fa ax 

While  the  sampling  rate  is  high  enough  for  the  complex  envelope 
representation,  the  signal  r(f)  usually  requires  a  higher  rate.  In  order  to  recover  the 
signal  from  the  Fourier  coefficients,  the  series  {AjJ  must  be  padded  with  zeros,  so 
that  the  modulated  signal  samples  occur  more  than  twice  per  carrier  period.  Both 


the  subroutines  ccq  and  calyce,  to  be  described  next,  use  six  samples  per  carrier 
period  for  plotting  purposes,  for  a  total  of  6 fcT0  samples  on  the  interval 
(-T0/2  ,  T0/ 2).  This  zero  padding  is  accomplished  by  using  the  value  Lc  =  6 fcT0 
instead  of  L  in  Eq.  (116). 

b.  Subroutine  Calyce 

The  subroutine  calyce  computes  the  output  electrical  signal,  both  in 
complex  envelope  y(t)  and  in  band  limited  y(t)  forms. 

(1)  Sampled  Output  Signal.  The  electrical  signal  ymn(t )  at  the 
output  of  the  element  (m,n)  in  the  receiver  array  is  given  by  Eq.  (40c). 
Substitution  of  the  expression  for  the  input  signal  transform,  Eq.  (118),  into  Eq. 
(40c)  yields 


w(k)  bk  H(fc  +  kfoJmn)  e^2*#7’*  l  ,  (119) 


k  =  -/„ 


ax 


after  performing  the  integration  indicated  in  Eq.  (40c),  and  where  f0T$  =  l/L  . 
From  the  definition  of  the  complex  envelope,  Eq.  (101),  the  sampled  complex 
envelope  of  the  output  signal  is  given,  by  inspection  of  Eq.  (119),  by 

^max 

V <m«T,)=  2)  «(t)  6k  Ha  +  .  (120) 

k  -  ax 


Equations  (119)  and  (120)  are  both  implemented  in  the  subroutine  calyce  The 
observations  made  about  the  recovery  of  the  signal  from  its  Fourier  coefficients — 
regarding  the  zero  padding  of  Fourier  coefficients — are  also  valid  for  y(t).  In  order 


to  have  a  correct  reconstruction  of  y(t),  the  constant  L  in  Eq.  (119)  is  substituted 
by  Lc  =  6j(.  T0,  corresponding  to  six  samples  per  carrier  period. 

(2)  Subroutine  Description.  Figure  9  lists  the  pseudo-code  of  the 
subroutine  calyce.  The  main  input  data  for  this  subroutine,  besides  those  listed  for 
the  subroutine  ccq,  are  a  three-dimensional  array  with  the  values  of  the  overall 
transfer  function  for  each  frequency  ( fc  +  kfj  and  position  (m,n)  in  the  receiver 
array,  and  the  Fourier  coefficients  bk  computed  in  the  subroutine  ccq. 

The  computation  of  the  toted  number  of  samples  L  and  sampling 
period  Ts  are  as  described  in  the  previous  section,  in  connection  to  the  generation  of 
the  transmitted  signal.  The  complex  envelope  samples  y^lJTs)  are  computed 
according  to  Eq.  (120).  The  signal  samples  y^lTs)  are  computed  using  Eq.  (119), 
with  zero  padding. 


Subroutine  Calyce 
Compute : 

total  number  of  samples  L  for  the  complex  envelope 
sample  period  TS  for  the  complex  envelope 
Compute  complex  envelope  samples  /Eq.  (120)/ 

Compute  and  plot  signal  /Eq.  (119)/ 

End  Subroutine  Calyce 


Fig.  9.  Simplified  Pseudo-Code  for  the  Subroutine  Calyce. 
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2.  Computing  the  Overall  Transfer  Function 

The  computation  of  the  overall  transfer  function  is  performed  by  a  group  of 
subprograms  subordinated  to  the  subroutine  chfmn ,  called  by  the  main  program. 
The  relationships  between  the  subroutines  used  for  this  implementation  are  shown 
in  Fig.  10.  The  program  was  implemented  in  FORTRAN,  with  most  of  the 
variables  shared  through  common  blocks . 

The  pseudo-code  for  the  subroutine  chfmn  is  shown  in  Fig.  11.  This 
subroutine  implements  Eq.  (100),  with  the  path  of  integration  F  truncated  to  the 
intervals  (-  jlower  ,  j  0)  plus  (0  ,  upper). 

3.  Integrand  Implementation 

With  reference  to  Fig.  10,  the  subprograms  related  to  the  computation  of 
the  integrand  are  reinteg,  iminteg,  integr,  refl,  hwvgky  and  dbsjO  *.  They  are  all 
implemented  as  FORTRAN  external  functions. 

The  functions  reinteg  and  iminteg  just  return  the  real  and  imaginary  parts 
of  the  integrand,  respectively.  They  are  called  by  the  integration  subroutine,  whose 
first  argument  is  the  name  of  the  function  to  be  integrated. 

a.  Function  Integr 

The  function  integr  computes  the  (complex)  integrand  in  Eq.  (100). 
The  pseudo-code  for  this  function  is  shown  in  Fig.  12.  The  argument  kv2  is 
interpreted  as  imaginary  when  it  has  a  negative  value,  in  all  subprograms.  The  IF 
statements  marked  in  the  list  are  related  to  the  medium  transfer  function.  New 


' 

'IMSL,  Inc.,  DBSJO ,  SFUN/  Library,  Bessel  Functions  of  Integer  Order,  1984. 
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Subroutine  Chfmn 

FOR  (  each  frequency  /=  X  +  k  jk  )  DO  /-A'max  Kmi.J 

EXECUTE  Weight  /compute  transmitter  array  weights  c^C /)/ 

FOR  (  each  position  yn  in  receiver  )  DO 

EXECUTE  Break  /compute  breakpoints  for  integration/ 
Compute  upper  and  lower  limits  of  integration 
FOR  (  each  position  xm  in  receiver  )  DO 

EXECUTE  QintegCReal  part)  /compute  Re{  H(j, rmn)  } 
and  its  error  Re{  EH(f,  rmn)  }/ 

EXECUTE  QintegCImag  part)  /compute  lm{  H(f, rmn)  } 
and  its  error  lm{  EH(f,  rmn)  }/ 

Compute  magnitude  and  phase  of  H  and  EH 
END  FOR  each  r_ 

m 

END  FOR  each  yn 

Plot  magnitude  and  phase  of  H  versus  position  (r,  y) 

END  FOR  each  frequency 
End  Subroutine  Chfnm 


Fig.  11.  Pseudo-Code  for  the  Subroutine  Chfmn. 
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Function  Integr(£y2) 

/this  function  is  executed  for  constant  frequency,  receiver 
position  ( m,n )  and  ky2/ 

Compute  kr  /Eq.  (92)/ 

IF1  (layered  waveguide)  THEN 

/?2,  =  Ref  1  (surface  data)  and  /?23  =  Refl  (bottom  data) 

END  IF 

/compute  product  'I'-£y2,  for  each  yq,  as  an  array  Ftraniq )/ 

FOR  (  each  position  yq  in  transmitter  array  )  DO 

IF1  (layered  waveguide)  Ftran(q)  =  Hwvgky  /  ky2  •  ^{f,ky2 ,  yQ)  yn)/ 
END  FOR  each  yq 

/compute  integrand  according  to  Eq.  (100)/ 

FOR  (each  position  xp  in  transmitter  array)  DO 

Compute  Bessel  function  J0(krR0)  /  R0  is  the  horizontal 

distance  from  transmitter  element  to  receiver  position/ 
FOR  (  each  position  yq  in  transmitter  array  )  DO 

Acumulate  product  cpq(/)  •  Ftraniq)  ■  J0(krR0)  /  S  X  / 
END  FOR  each  y 
END  FOR  each  xp 

Integr  =  (acumulated  product)/2r 
End  Function  Integr 

■Other  IF  statements  may  be  necessary  for  different  medium  transfer 
functions 

Fig.  1  i.  Pseudo-Code  for  the  Function  Integr. 


transfer  functions  are  implemented  by  adding  new  IF  statements.  The  functions 
integr  and  break,  to  be  described  later,  are  the  only  ones  that  need  to  be  modified 
when  incorporating  other  medium  transfer  functions  to  the  program. 

b.  Function  Refl 

This  function  computes  the  complex  reflection  coefficients.  Its 
pseudo -code  is  shown  in  Fig.  13. 

The  expressions  for  both  reflection  coefficients,  Eqs.  (73)  and  (76),  can 
be  written  as 


Pp  ^y2 

Pp  ky2 


kyp 
Pi  kyp 


(121) 


Function  ftefl(cp,  pp,  Ay2) 

/computes  reflection  coefficient  Rip/ 

Compute  kyp  /Eq.  (92)/ 

IF  (the  speeds  of  sound  are  different)  THEN 

Compute  reflection  coefficient  by  Eq. (121) 
ELSE 

Compute  reflection  coefficient  by  Eq.  (122) 
END  IF 

End  Function  Refl 


Fig.  13.  Pseudo-Code  for  the  Function  Refl. 
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where  kyp  is  given  by  Eq.  (92).  If  the  speed  of  sound  is  the  same  in  both  media, 
then  Eq.  (121)  can  be  simplified  to 


Pp  ~  Pi 
Pp  +  Pi 


(122) 


c.  Function  Hwvgky 

The  function  hwvgky  returns  the  complex  product  ky2  ^ {lKz,y0‘,ya)  , 
the  medium  dependent  factor  in  the  integral  expression  for  the  overall  transfer 
function,  Eq.  (100).  This  function  is  specific  to  the  layered  waveguide  and 
implements  Eq.  (81),  multiplied  by  iy2.  It  is  called  by  the  function  integr,  as 
described  earlier  (see  Figs.  10  and  12).  Its  implementation  is  trivial.  Its  arguments 
axe  the  frequency,  ky2  and  the  y  coordinates  of  the  transmitter  and  receiver  array 
elements.  y  and  yn  ,  respectively.  The  reflection  coefficients  /£>,  and  /?23,  and  the 
depth  D  are  available  through  common  blocks. 

d.  Examples 

In  this  section  we  examine  the  behavior  of  the  integrand  in  the 
expression  for  the  overall  transfer  function,  Eq.  (100).  In  all  plots  in  which  ky.2  is 
the  independent  variable,  negative  values  of  ky.2  are  imaginary  and  positive  values 
are  real.  For  simplicity,  both  negative  imaginary  and  positive  real  values  of  k  „  are 
plotted  on  the  same  axis.  The  symbol  appearing  on  the  plots  is  used  by  the 
plotting  routine  to  mark  different  plots  on  the  same  graph,  but  it  has  no  special 
meaning  in  the  present  context. 

(1)  The  Flange  Factor  J0{krR0).  Figures  14  and  15  show  the  plots  of 
J0(krRo)  vs  £y2,  which  is  related  to  kr  through  Eq.  (91).  The  frequency  is  /= 
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Fig.  16.  Range  Factor  JJkrR0)  for  /=  600  Hz,  c,  =  1600  m/s,  and 

R0  =  1000  m  (cont.). 
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=  600  Hz  and  the  sound  speed  is  Cj  =  1,500  m/s,  which  results  in  a  wavenumber 
*2  =  2.09  m~*  .  In  Fig.  14,  the  distance  is  R0  =  100  m  and,  in  Fig.  15,  it  is  1,000  m. 
In  those  figures  the  range  of  ky2  is  (-j’0.25  ,  jO  )  plus  (  0  ,  ^  .  An  important 

characteristic  of  the  factor  J0(kTR0)  is  that  it  oscillates  slower  near  ky2  a  0  . 
Another  characteristic  behavior  is  that  it  oscillates  faster  as  the  range  increases.  If 
the  other  terms  in  the  integrand  also  have  "slow"  oscillations,  then  the  greater 
contribution  for  the  result  of  the  integration  comes  from  the  region  ky2  a  0  ,  as  a 
result  of  the  smoothing  nature  of  the  integration  operation.  Physically,  the  plane 
waves  that  travel  near  the  horizontal,  that  is,  those  with  ky2  ft  0  ,  contribute  more 
to  the  total  field  than  those  that  interact  more  with  the  boundaries.  In  order  to 
compare  the  "frequency  of  oscillations"  of  the  Bessel  function  and  of  the  medium 
transfer  function,  recall  that,  for  krRlf  >  2jt.  the  zeros  of  J0(krR0)  occur 
approximately  at  intervals  of 

6k=  —  .  (123) 

Ro 

(2)  The  Medium  Factor.  Figure  16  shows  the  plot  of  the  term  ky2 
for  the  layered  waveguide  whose  characteristics  are  shown  in  Table  1.  The  depth  D 
of  the  ocean  is  100  m  and  the  frequency  is  500  Hz.  The  position  of  thi  point  source 
is  ys ,  ^s)  =  (0,  30,  0)  meters  and  the  point  receiver  sensor  is  at  (rT  ,yr  .  cr)  = 
=  (1000  ,50  ,  0)  meters,  corresponding  to  the  range  factor  shown  in  Fig.  15. 

This  waveguide  is  of  the  almost  -pressure -release -surface/fast  - 
bottom  type.  The  peaks  are  evident  in  the  figures.  From  the  discussion  in  Section 
III -A -lb,  there  are  33  peaks  in  the  range  0  <  A:y2  <  1  -04  m  ^  ,  separated  by 
Afcy2  a  it/  100. 
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TABLE  1.  WAVEGUIDE  ACOUSTIC  CHARACTERISTICS 


Medium 

Speed  of  Sound 

Density 

(m  s-1) 

(Kg  nO3 

I 

343.0 

1.21 

II 

1500.0 

1026.0 

III 

1730.0 

2070.0 

A  criterion  to  set  an  upper  limit  of  integration  in  Eq.  (100) — the 
upper  discussed  in  Section  III-B-2 — could  be  based  on  the  comparison  between  the 
"periodicities"  of  the  medium  and  range  factors  in  the  integrand.  If  the  factors  in 
the  expression  fur  Eq.  (81),  are  multiplied  out,  the  faster  oscillatory  factor — the 
one  with  the  greater  exponent — has  a  "period"  of  x/Dor  more,  with  respect  to  £y2. 
that  is, 

Sk  z  —  .  (124) 

D 

From  Eq.  (91),  the  increments  Sky2  and  Skr  are  approximately  related  by 


when  ky2  is  real.  Equations  (123)  and  (125)  can  be  used  to  compute  the  interval 
between  zeros  of  the  range  factor,  when  it  is  plotted  versus  ky2.  If  we  require,  for 
example,  that  the  integral  be  truncated  when  the  Bessel  function  passes  through 
zero  N times  during  one  "period"  of  the  medium  factor,  then,  substituting  Eq.  (123) 
and  Sky2  =  r  /  N  D  into  Eq.  (125),  yields 

k  2 

uppers—  - - - —  —  .  (126) 

Jl  +  ( Ro/ND )2 

The  upper  limit  computed  according  to  Eq.  (126)  can  smaller  than  the  last  peak 

f—n - j 

position,  given  by  v  ^2  _  ^3  [see  discussion  following  Eq.  (99)].  In  this  case,  the 
limit  could  be  increased,  or  those  peaks  outside  the  limit  could  be  simply  discarded. 
In  this  implementation,  the  upper  limit  is  set  to  L,. 

The  lower  limit  of  integration,  in  the  region  where  ky2  is 
imaginary  negative  and  |  ty2-v[r  |  decreases  exponentially  as  k  2-  -jtc  .  should  be 
set  to  minimize  the  error  of  integration  for  the  slower  decaying  factor  in  &v2-' V — the 
one  with  the  exponent  of  smallest  absolute  value.  The  worst  case  is  when  there  is 
one  constant  term,  which  happens  when  the  source  and  the  receiver  are  at  the  same 
depth,  for  example.  Note  that  the  term  ky2-^>  is  pure  imaginary  when  ky2  is 
negative  imaginary,  and  decreases  towards  zero  as  e  I^y2i J/o  _  f  ~0|£y2|  The  lower 
limit  is  set  to  -j0.125  k^  in  this  implementation. 

(3)  The  Integrand.  Figure  17  shows  the  plot  of  the  integrand  for  the 
waveguide  discussed  above,  at  /=  500  Hz,  for  a  range  =  1,000  m.  For  ky2 
greater  than  about  1.5,  corresponding  to  Ns  10  in  Eq.  (126),  the  integrand  is  fairly 
symmetrical  around  zero,  suggesting  a  low  value  for  the  integral  in  that  region. 


76 


4.  Other  Subroutines 
a.  Subroutine  Break 

This  subroutine  divides  the  interval  of  integration  in  order  to  decrease 
the  error  of  the  integration  subroutine.  The  criteria  used  in  this  division  are  the 
behavior  of  the  medium  transfer  function  and  the  relationship  between  the 
"periodicities"  of  that  function  and  of  the  Bessel  function  in  the  integrand.  For 
each  different  medium  transfer  function,  a  different  interval  subdivision  strategy 
must  be  implemented.  The  upper  and  lower  limits  of  integration  are  independent  of 
the  type  of  medium,  in  this  implementation.  Nevertheless,  they  could  be  made 
medium  dependent,  in  which  case  their  computation  should  be  included  in  the 
subroutine  break.  Figure  18  lists  the  pseudo-code  for  this  subroutine. 

In  the  case  of  the  lavered  waveguide  function,  the  breakpoints  are  set 
equal  to  the  locations  of  the  peaks,  if  any,  as  discussed  in  Section  III -A -lb, 
regarding  the  almost  -  pressure -release -surface/fast -bottom.  The  interval  of  values 
of  ky2  outside  the  range  of  the  peaks  is  divided  evenly  into  subintervals  of  size 
b*  /  D ,  corresponding  to  five  "periods,"  as  given  by  Eq.  (124). 

(1)  Subroutine  Description.  The  input  parameters  for  the  subroutine 
break  are  the  frequency  J.  the  speed  of  sound  a, — computed  in  the  subroutine  chjmn 
and  taking  into  account  the  gradient  of  the  speed  of  sound,  for  the  case  of  other 
media.  Other  variables  shared  thiough  common  blocks  are  the  type  of  medium. 
OCNTF ,  and  the  waveguide  characteristics.  The  outputs  are  the  number  of 
breakpoints  and  a  vector  with  the  values  of  kyi  at  the  breakpoints. 

Initially,  the  subroutine  divides  the  interval  of  integration  evenly. 
Then,  it  proceeds  to  check  for  the  presence  of  peaks  in  the  medium  transfer 
function  as  discussed  in  Section  III -A -lb.  The  peaks  are  the  zeros  of  Eq.  (99). 
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Subroutine  Break(/,  Cj  ,  breakpoints  ,  number_brkpts) 

IF<  (layered  waveguide)  THEN 

/divide  range  evenly  (-0.125  k 2,  0)  U  (0,  k^)/ 
number Jbrkpts  -  (1.125  k^  r  A kyi)  -  1 
FOE  (n  =  1  to  number_brkpts)  DO 

breakpoints^ n)  =  (—0.125  k^)  +  n  Afcy2 

IF  ( breakpointsin )  <  0)  index  =  n  /index  keeps  track  of 
last  nonpositive  element  of  vector  breakpoints/ 

END  FOR  n 

/check  for  presence  of  peaks/ 

IF  (peaks  do  exist)  THEN  /a>  xf2  ,  see  Eq.  (99)/ 

/the  peaks  are  located  at  &y2  =  0n ±D,  see  Eq.  (99)/ 
number_peaks  =  [(  a/x)  -  0 . 5]  /  discussion  after  Eq.  (99)/ 

FOR  (n=  indei+l  to  tnder+l+  number  peaks )  DO 
Compute  n_th_zero  of  the  function  Denwgd 
breakpoint sin)  =  n_th_:ero  ?  D 
END  FOR  n 

Update  number_brkpts  to  reflect  total  number  of  breakpoints 
END  IF 
END  IF 

End  Subroutine  Break 

•Other  IF  statements  may  be  necessary  for  different  medium  transfer 
functions 

Fig.  lb.  Pseudo-Code  for  the  Subroutine  Break 


implemented  as  the  FORTRAN  external  function  denwgd.  The  zeros  are  found  by 
the  subroutine  dzbren1.  The  (positive)  breakpoints  originally  computed  are  changed 
to  the  values  of  the  positions  of  the  peaks.  This  implementation  assumes  a  lowei 
limit  of  integration  of  -j  0.125  k%  and  an  upper  limit  equal  to  k 2. 
b.  Subroutine  Weight 2 

The  subroutine  weight  computes  the  transmitter  array  complex  weights 
Cpq(/).  The  input  parameters  are  the  frequency;  the  transmitter  array  data;  the 
direction  cosines  u'  (with  respect  to  the  X  direction)  and  v1  (with  respect  to  the  Y 
direction)  of  the  transmitter  array  main  lobe;  the  speed  of  sound  at  the  center  of  the 
array — set  to  c-%  for  the  layered  waveguide,  but  included  in  the  present 
implementation  to  allow  for  other  media — and  the  gradient  of  the  speed  of 
sound — set  to  zero  in  this  implementation.  The  complex  weights  are  given  by  [Ref. 
5:  p.  133] 


£pq  if)  —  apq 


(127) 


where 


flpq  —  Ap  Bq  , 


(126) 


7pq 


( 129) 


'IMSL,  Inc.,  DZBREN,  MA  TH/ Library,  Nonlinear  Equations ,  1985. 

cThis  implementation  is  based  on  a  program  originally  written  by  Dr. 
Lawrence  J.  Ziomek.  at  the  Naval  Postgraduate  School. 

c>3 


(130) 


a  =  5K) 

/ 

Ap  and  Bq  axe  the  amplitude  weights  along  the  X  and  Y  directions,  respectively,  X 
is  the  (depth-dependent)  wavelength,  dxt  and  dyt  are  the  interelement  spacing 
along  the  X  and  Y  directions,  respectively,  and  c^iy^)  is  the  speed  of  sound  at  the 
depth  yq. 

c.  Subroutine  Qinteg 

The  subroutine  qinteg  is  an  interface  for  the  subroutine  dqSags *,  which 
actually  performs  the  integration  in  Eq.  (100).  The  input  arguments  are  the  name 
of  the  external  FORTRAN  function  to  be  integrated,  the  limits  of  integration,  the 
breakpoints,  the  number  of  breakpoints,  and  the  required  error  in  the  result.  The 
output  arguments  are  the  result  of  the  integration  and  its  estimated  error. 

Initially,  each  interval  of  integration,  defined  by  the  breakpoints,  is 
subdivided  in  order  to  improve  the  convergence  of  the  subroutine  dq2ags ,  which  is 
called  by  qinteg  to  perform  each  intermediate  integration.  If  the  error  of  the 
integration  of  a  subinterval,  as  computed  by  dq2ags,  exceeds  the  requirements,  that 

l 

subinterval  is  furlhei  subdivided.  Finally,  the  results  and  errors  of  integration 
provided  by  dqSags  are  accumulated. 

As  indicated  in  the  pseudo-code  of  the  subroutine  chfmn,  Fig.  11.  each 


•IMSL.  Inc.,  DQ2AGS ,  MATH/ Library.  Integration  and  Differentiation , 

1985. 
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integration  is  performed  twice,  for  the  real  and  imaginary  parts  of  the  integrand.  In 
terms  of  computing  time,  this  is  the  most  critical  section  of  code.  Not  only  does  the 
subroutine  qinteg  calls  dqSags  many  times — at  least  ten  times  between  breakpoints, 
in  the  present  implementation — but,  what  is  worse,  dqSags  c an  execute  the  function 
to  be  integrated,  ultimately  integr ,  hundreds  of  times.  The  computation  .ime 
increases  with  frequency  and  range: 

•  The  interval  of  integration  increases  as  0(1 i^),  directly  proportional  to 
frequency. 

•  As  the  range  increases,  the  range  factor  J0(krR0)  oscillates  faster,  degrading  the 
convergence  of  the  integration  subroutine,  causing  more  executions  of  the 
function  integr  and  more  frequent  interv  .1  subdivisions. 

The  effect  of  range  can  be  minimized  if  the  limits  of  integration,  lower  and  upper, 
decrease  (in  absolute  value)  with  range,  as  in  Eq.  (126)  for  the  upper  limit. 

C.  RESULTS 

The  algorithms  described  in  the  previous  section  are  tested  by  using  two 
examples  for  which  the  results  are  easily  interpreted.  The  unbounded  homogeneous 
medium  was  presented  in  Section  II-C-3a  as  a  particular  case  of  the  layered 
waveguide,  all  three  media  with  the  same  acoustic  characteristics.  The  surface 
reflection  problem  was  analyzed  in  Section  II -C -3b.  In  both  cases,  the  input  to  the 
system  is  a  time-harmonic  point  source  and  the  output  electrical  signal — 
numerically  equal  to  the  velocity  potential — is  given  by  Eqs.  (85)  and  (88).  From 
Ec  40c).  the  time-independent  term  of  the  output  signal  is,  for  a  time-harmonic 
input,  equal  to  the  overall  transfer  function,  that  is, 
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(40c) 


*„(0  =  j“xu)  j'2’! 1  df , 

-oo 

!U(0  =  ["«(/-/„)  WUr„)  ^'2T/1  df,  (131a) 

-oo 

«,„(<)  -HarJ*’’2'4'.  (131b) 

Therefore,  the  magnitude  of  the  overall  transfer  function  is  equal  to  the  magnitude 
of  the  velocity  potential. 

Results  for  the  waveguide  example  used  so  far  (see  Table  1)  are  also  presented. 

i  ii  /-»*-►-»  r\n-on  strut  «* 

a.  iiuiiiu^cucutio  lViculUill 

Figure  19  presents  the  plot  of  the  magnitude  of  the  overall  transfer  function 
as  a  function  of  horizontal  distance  r  from  the  source.  The  acoustic  characteristics 
are  those  shown  in  Table  1  for  medium  II.  The  depth  of  the  source  is  30  m  and  of 
the  receiver  array,  50  m.  The  expected  result  is  a  straight  line  on  a  log-log  plot, 
corresponding  to  a  velocity  potential  proportional  to  the  inverse  of  distance.  At  100 
m,  the  error  of  the  computed  value  is  6%,  and  increases  with  distance  up  to  J9LA  at 
3,100  m.  The  integration  routine  computed  the  integrand  3,360  times  (total  for 
real  and  imaginary  parts)  for  100  m  of  distance,  a  number  that,  increased  to  34.776 
at  3,100  m,  indicating  an  increased  difficulty  in  convergence  of  the  integral.  An 
analogous  computation  for  the  layered  waveguide  took  45,000  and  70.000 
computations  of  the  integrand,  for  100  n.  and  3.100  m,  respectively. 
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2.  Surface  Reflection 

Figure  20  shows  the  overall  transfer  function  as  a  function  of  depth,  for  a 
distance  of  1,000  m  between  source  and  receiver.  The  acoustic  characteristics  of  the 
two  media  are  those  shown  in  Table  1  for  media  I  and  II.  The  surface  is,  for  all 
practical  purposes,  pressure  release.  The  depth  of  the  source  is  30  m.  For  these 
data,  the  approximations  for  Eq.  (89)  are  valid,  indicating  a 
with  depth.  The  error  increases  from  zero  at  the  surface,  up  to 
30  m,  the  depth  of  the  source — see  discussion  in  Section  III -B -3d,  regarding  the 
lower  Lnit  of  integratic.i.  At  30  m,  near  a  theoretical  maximum,  the  error  is  45%  . 
At  depths  greater  than  60  m,  the  error  is  negligible.  For  comparison  purposes,  the 
integrand  was  computed  about  11,000  times  for  each  value  of  the  overall  transfer 
function. 

3.  Layered  Waveguide:  Waveform  Prediction 

Figure  21  shows  the  transmitted  signal.  This  signal  was  obtained  by 
truncating  the  Fourier  series  of  a  CW  pulse,  as  explained  in  Section  III -B -la.  The 
reference  signal  is  a  rectangular -envelope  CW  pulse  of  duration  20  ms  and 
repetition  period  of  400  ms.  The  carrier  frequency  is  jf  =  500  Hz  Seventeen 
frequency  components  are  used  for  the  transmitted  signal,  that  is,  Kmav  =  8.  The 
resultant  bell -shaped  CW  pulse  has  a  duration  of  100  ms.  Ringing  is  barely  visible, 
due  to  the  low  sidelobes  of  the  Hamming  window  used  for  the  truncation  of  the 
Fourier  series.  The  extension  of  time  in  the  plot  of  Fig.  21  is  equal  to  one  period. 


siryjsoidal  variation 


a  ijiaximum  error  at 
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Fig.  20.  Overall  Transfer  Function  for  the  Surface  Reflection 
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The  acoustic  characteristics  of  the  media  are  the  same  as  for  the  examples 
in  Section  III -B -3d  and  are  listed  in  Table  1.  The  distance  is  R0  =  1,000  m.  The 
transmitter  is  located  at  (j^,  ys ,  2s)  =  (0,  30,  0)  meters,  and  the  receiver  at 
(xT,  yT ,  2r)  =  (1000,  50,  0)  meters. 

Figure  22  shows  the  received  pulse  when  the  transmitter  is  a  single  point 
source.  The  distortion  due  to  the  multipath  nature  of  the  medium  is  evident. 

Figure  23  shows  the  received  pulse  due  to  a  5*5  transmitter  array.  The 
interelement  spacing  is  roughly  A/2  in  both  the  X  and  Y  directions.  Note  that  the 
pulse  arrivals  occur  at  the  same  instants  as  in  Fig.  22,  but  the  relative  amplitudes 
are  different. 

D.  CONCLUSIONS 

The  equations  that  describe  a  range -independent  acoustic  communication 
channel  were  derived  by  using  linear  systems  theory,  a  basic  engineering  tool,  as  a 
framework.  They  incorporate  aspects  of  signal  processing  and  acoustic  propagation. 
The  main  results  from  Chapter  II  are  the  expression  for  the  overall  transfer 
function,  Eq.  (51)  or  its  alternative  form.  Eq.  (60):  and  the  procedure  to  obtain  the 
transformed  medium  transfer  function,  the  solution  to  the  inhomogeneous 
Helmholtz  wave  equation.  Eq.  (90).  These  results,  together  with  Eq.  (40c).  were 
used  to  implement  a  waveform  prediction  program,  as  described  in  Section  III -B. 
We  do  not  claim  novelty  of  results.  Equation  (60)  is  a  mere  extension  of  a  classical 
result,  in  which  the  contributions  of  the  point  sources  that  make  up  an  array  are 


Fig.  22.  Output  Electrical  Signal  Due  to  a  Single  Point  Source. 
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simply  added,  a  consequence  of  the  linearity  of  the  wave  equation  that  describes 
small -amplitude  acoustic  phenomena.  The  medium  transfer  function  for  the  layered 
waveguide  was  derived  in  order  to  enable  us  to  implement  a  working  computer 
program. 

The  program  described  in  Section  III -B  is  based  on  a  slightly  modified  form  of 
Eq.  (60),  in  which  we  use  the  wavenumber  y-component  iy2  as  the  independent 
variable.  The  only  justification  for  this  change  of  variable  is  the  easier  analysis  of 
the  medium  transfer  function.  The  program  works,  but  is  slow.  The  computation 
of  the  waveform  shown  in  Fig.  23,  for  example,  took  about  three  minutes  per 
frequency  component,  with  an  overhead  time — independent  of  the  number  of 
frequency  components — of  about  two  minutes,  in  an  IBM  3033U/4381  2-cpu 
network.  Both  the  speed  and  accuracy  can  be  improved  by  using  an  adaptive 
technique  for  the  truncation  of  the  interval  of  integration,  that  is,  computation  of 
the  limits  of  integration  as  a  function  of  the  behavior  of  the  medium  and  range 
factors  in  the  integrand. 

Other  implementations  are  possible,  for  example,  the  Fast -Field -Program 
described  in  [Ref.  4:pp.  90-92].  Equation  (51).  in  the  form  of  a  two-dimensional 
Fourier  transform,  could  be  useful  for  implementation. 

Concerning  the  medium  transfer  function,  analytical  solutions  to  the 
range -independent  wave  equation  are  only  available  for  a  few  speed -of -sound 
profiles.  Approximate  solutions,  like  those  provided  by  the  WKB  method,  are 
candidates  for  implementation.  On  the  other  hand,  numerical  methods  may  not  be 
feasible  for  our  purpose,  waveform  prediction,  due  to  time  constraints. 
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As  mentioned  in  the  beginning  of  Chapter  III,  the  main  thrust  for  the  present 
implementation  of  the  range -independent  equations  was  to  have  a  working  program. 
The  program  is  working.  More  than  an  improvement,  it  possibly  needs  a  different 
implementation.  Nevertheless,  the  basic  important  fact  is  that,  independent  of  the 
implementation,  the  modular  nature  of  the  equations  should  be  fully  exploited. 
Solutions  to  specific  problems  can  be  implemented  by  adding  new  subprograms, 
with  a  minimum  of  code  change.  Further  investigation  is  recommended. 
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